home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / meschach / meschach1.shar < prev    next >
Encoding:
Text File  |  1994-06-07  |  139.1 KB  |  5,841 lines  |  [TEXT/ttxt]
  1. # to unbundle, sh this file (in an empty directory)
  2. echo copy.c 1>&2
  3. sed >copy.c <<'//GO.SYSIN DD copy.c' 's/^-//'
  4. -
  5. -/**************************************************************************
  6. -**
  7. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  8. -**
  9. -**                 Meschach Library
  10. -** 
  11. -** This Meschach Library is provided "as is" without any express 
  12. -** or implied warranty of any kind with respect to this software. 
  13. -** In particular the authors shall not be liable for any direct, 
  14. -** indirect, special, incidental or consequential damages arising 
  15. -** in any way from use of the software.
  16. -** 
  17. -** Everyone is granted permission to copy, modify and redistribute this
  18. -** Meschach Library, provided:
  19. -**  1.  All copies contain this copyright notice.
  20. -**  2.  All modified copies shall carry a notice stating who
  21. -**      made the last modification and the date of such modification.
  22. -**  3.  No charge is made for this software or works derived from it.  
  23. -**      This clause shall not be construed as constraining other software
  24. -**      distributed on the same medium as this software, nor is a
  25. -**      distribution fee considered a charge.
  26. -**
  27. -***************************************************************************/
  28. -
  29. -
  30. -static    char    rcsid[] = "$Id: copy.c,v 1.2 1994/01/13 05:37:14 des Exp $";
  31. -#include    <stdio.h>
  32. -#include    "matrix.h"
  33. -
  34. -
  35. -
  36. -/* _m_copy -- copies matrix into new area */
  37. -MAT    *_m_copy(in,out,i0,j0)
  38. -MAT    *in,*out;
  39. -u_int    i0,j0;
  40. -{
  41. -    u_int    i /* ,j */;
  42. -
  43. -    if ( in==MNULL )
  44. -        error(E_NULL,"_m_copy");
  45. -    if ( in==out )
  46. -        return (out);
  47. -    if ( out==MNULL || out->m < in->m || out->n < in->n )
  48. -        out = m_resize(out,in->m,in->n);
  49. -
  50. -    for ( i=i0; i < in->m; i++ )
  51. -        MEM_COPY(&(in->me[i][j0]),&(out->me[i][j0]),
  52. -                (in->n - j0)*sizeof(Real));
  53. -        /* for ( j=j0; j < in->n; j++ )
  54. -            out->me[i][j] = in->me[i][j]; */
  55. -
  56. -    return (out);
  57. -}
  58. -
  59. -/* _v_copy -- copies vector into new area */
  60. -VEC    *_v_copy(in,out,i0)
  61. -VEC    *in,*out;
  62. -u_int    i0;
  63. -{
  64. -    /* u_int    i,j; */
  65. -
  66. -    if ( in==VNULL )
  67. -        error(E_NULL,"_v_copy");
  68. -    if ( in==out )
  69. -        return (out);
  70. -    if ( out==VNULL || out->dim < in->dim )
  71. -        out = v_resize(out,in->dim);
  72. -
  73. -    MEM_COPY(&(in->ve[i0]),&(out->ve[i0]),(in->dim - i0)*sizeof(Real));
  74. -    /* for ( i=i0; i < in->dim; i++ )
  75. -        out->ve[i] = in->ve[i]; */
  76. -
  77. -    return (out);
  78. -}
  79. -
  80. -/* px_copy -- copies permutation 'in' to 'out' */
  81. -PERM    *px_copy(in,out)
  82. -PERM    *in,*out;
  83. -{
  84. -    /* int    i; */
  85. -
  86. -    if ( in == PNULL )
  87. -        error(E_NULL,"px_copy");
  88. -    if ( in == out )
  89. -        return out;
  90. -    if ( out == PNULL || out->size != in->size )
  91. -        out = px_resize(out,in->size);
  92. -
  93. -    MEM_COPY(in->pe,out->pe,in->size*sizeof(u_int));
  94. -    /* for ( i = 0; i < in->size; i++ )
  95. -        out->pe[i] = in->pe[i]; */
  96. -
  97. -    return out;
  98. -}
  99. -
  100. -/*
  101. -    The .._move() routines are for moving blocks of memory around
  102. -    within Meschach data structures and for re-arranging matrices,
  103. -    vectors etc.
  104. -*/
  105. -
  106. -/* m_move -- copies selected pieces of a matrix
  107. -    -- moves the m0 x n0 submatrix with top-left cor-ordinates (i0,j0)
  108. -       to the corresponding submatrix of out with top-left co-ordinates
  109. -       (i1,j1)
  110. -    -- out is resized (& created) if necessary */
  111. -MAT    *m_move(in,i0,j0,m0,n0,out,i1,j1)
  112. -MAT    *in, *out;
  113. -int    i0, j0, m0, n0, i1, j1;
  114. -{
  115. -    int        i;
  116. -
  117. -    if ( ! in )
  118. -    error(E_NULL,"m_move");
  119. -    if ( i0 < 0 || j0 < 0 || i1 < 0 || j1 < 0 || m0 < 0 || n0 < 0 ||
  120. -     i0+m0 > in->m || j0+n0 > in->n )
  121. -    error(E_BOUNDS,"m_move");
  122. -
  123. -    if ( ! out )
  124. -    out = m_resize(out,i1+m0,j1+n0);
  125. -    else if ( i1+m0 > out->m || j1+n0 > out->n )
  126. -    out = m_resize(out,max(out->m,i1+m0),max(out->n,j1+n0));
  127. -
  128. -    for ( i = 0; i < m0; i++ )
  129. -    MEM_COPY(&(in->me[i0+i][j0]),&(out->me[i1+i][j1]),
  130. -         n0*sizeof(Real));
  131. -
  132. -    return out;
  133. -}
  134. -
  135. -/* v_move -- copies selected pieces of a vector
  136. -    -- moves the length dim0 subvector with initial index i0
  137. -       to the corresponding subvector of out with initial index i1
  138. -    -- out is resized if necessary */
  139. -VEC    *v_move(in,i0,dim0,out,i1)
  140. -VEC    *in, *out;
  141. -int    i0, dim0, i1;
  142. -{
  143. -    if ( ! in )
  144. -    error(E_NULL,"v_move");
  145. -    if ( i0 < 0 || dim0 < 0 || i1 < 0 ||
  146. -     i0+dim0 > in->dim )
  147. -    error(E_BOUNDS,"v_move");
  148. -
  149. -    if ( (! out) || i1+dim0 > out->dim )
  150. -    out = v_resize(out,i1+dim0);
  151. -
  152. -    MEM_COPY(&(in->ve[i0]),&(out->ve[i1]),dim0*sizeof(Real));
  153. -
  154. -    return out;
  155. -}
  156. -
  157. -/* mv_move -- copies selected piece of matrix to a vector
  158. -    -- moves the m0 x n0 submatrix with top-left co-ordinate (i0,j0) to
  159. -       the subvector with initial index i1 (and length m0*n0)
  160. -    -- rows are copied contiguously
  161. -    -- out is resized if necessary */
  162. -VEC    *mv_move(in,i0,j0,m0,n0,out,i1)
  163. -MAT    *in;
  164. -VEC    *out;
  165. -int    i0, j0, m0, n0, i1;
  166. -{
  167. -    int        dim1, i;
  168. -
  169. -    if ( ! in )
  170. -    error(E_NULL,"mv_move");
  171. -    if ( i0 < 0 || j0 < 0 || m0 < 0 || n0 < 0 || i1 < 0 ||
  172. -     i0+m0 > in->m || j0+n0 > in->n )
  173. -    error(E_BOUNDS,"mv_move");
  174. -
  175. -    dim1 = m0*n0;
  176. -    if ( (! out) || i1+dim1 > out->dim )
  177. -    out = v_resize(out,i1+dim1);
  178. -
  179. -    for ( i = 0; i < m0; i++ )
  180. -    MEM_COPY(&(in->me[i0+i][j0]),&(out->ve[i1+i*n0]),n0*sizeof(Real));
  181. -
  182. -    return out;
  183. -}
  184. -
  185. -/* vm_move -- copies selected piece of vector to a matrix
  186. -    -- moves the subvector with initial index i0 and length m1*n1 to
  187. -       the m1 x n1 submatrix with top-left co-ordinate (i1,j1)
  188. -        -- copying is done by rows
  189. -    -- out is resized if necessary */
  190. -MAT    *vm_move(in,i0,out,i1,j1,m1,n1)
  191. -VEC    *in;
  192. -MAT    *out;
  193. -int    i0, i1, j1, m1, n1;
  194. -{
  195. -    int        dim0, i;
  196. -
  197. -    if ( ! in )
  198. -    error(E_NULL,"vm_move");
  199. -    if ( i0 < 0 || i1 < 0 || j1 < 0 || m1 < 0 || n1 < 0 ||
  200. -     i0+m1*n1 > in->dim )
  201. -    error(E_BOUNDS,"vm_move");
  202. -
  203. -    if ( ! out )
  204. -    out = m_resize(out,i1+m1,j1+n1);
  205. -    else
  206. -    out = m_resize(out,max(i1+m1,out->m),max(j1+n1,out->n));
  207. -
  208. -    dim0 = m1*n1;
  209. -    for ( i = 0; i < m1; i++ )
  210. -    MEM_COPY(&(in->ve[i0+i*n1]),&(out->me[i1+i][j1]),n1*sizeof(Real));
  211. -
  212. -    return out;
  213. -}
  214. //GO.SYSIN DD copy.c
  215. echo err.c 1>&2
  216. sed >err.c <<'//GO.SYSIN DD err.c' 's/^-//'
  217. -
  218. -/**************************************************************************
  219. -**
  220. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  221. -**
  222. -**                 Meschach Library
  223. -** 
  224. -** This Meschach Library is provided "as is" without any express 
  225. -** or implied warranty of any kind with respect to this software. 
  226. -** In particular the authors shall not be liable for any direct, 
  227. -** indirect, special, incidental or consequential damages arising 
  228. -** in any way from use of the software.
  229. -** 
  230. -** Everyone is granted permission to copy, modify and redistribute this
  231. -** Meschach Library, provided:
  232. -**  1.  All copies contain this copyright notice.
  233. -**  2.  All modified copies shall carry a notice stating who
  234. -**      made the last modification and the date of such modification.
  235. -**  3.  No charge is made for this software or works derived from it.  
  236. -**      This clause shall not be construed as constraining other software
  237. -**      distributed on the same medium as this software, nor is a
  238. -**      distribution fee considered a charge.
  239. -**
  240. -***************************************************************************/
  241. -
  242. -
  243. -/*
  244. -  File with basic error-handling operations
  245. -  Based on previous version on Zilog
  246. -  System 8000 setret() etc.
  247. -  Ported to Pyramid 9810 late 1987
  248. -  */
  249. -
  250. -static    char    rcsid[] = "$Id: err.c,v 1.5 1994/01/13 05:37:35 des Exp $";
  251. -
  252. -#include    <stdio.h>
  253. -#include    <setjmp.h>
  254. -#include    <ctype.h>
  255. -#include        "err.h"
  256. -
  257. -
  258. -#ifdef SYSV
  259. -/* AT&T System V */
  260. -#include    <sys/signal.h>
  261. -#else
  262. -/* something else -- assume BSD or ANSI C */
  263. -#include    <signal.h>
  264. -#endif
  265. -
  266. -
  267. -
  268. -#define        FALSE    0
  269. -#define        TRUE    1
  270. -
  271. -#define    EF_EXIT        0
  272. -#define    EF_ABORT    1
  273. -#define    EF_JUMP        2
  274. -#define    EF_SILENT    3
  275. -
  276. -/* The only error caught in this file! */
  277. -#define    E_SIGNAL    16
  278. -
  279. -static    char    *err_mesg[] =
  280. -{      "unknown error",            /* 0 */
  281. -      "sizes of objects don't match",    /* 1 */
  282. -      "index out of bounds",        /* 2 */
  283. -      "can't allocate memory",        /* 3 */
  284. -      "singular matrix",            /* 4 */
  285. -      "matrix not positive definite",    /* 5 */
  286. -      "incorrect format input",        /* 6 */
  287. -      "bad input file/device",        /* 7 */
  288. -      "NULL objects passed",        /* 8 */
  289. -      "matrix not square",            /* 9 */
  290. -      "object out of range",        /* 10 */
  291. -      "can't do operation in situ for non-square matrix",   /* 11 */
  292. -      "can't do operation in situ",        /* 12 */
  293. -      "excessive number of iterations",    /* 13 */
  294. -      "convergence criterion failed",    /* 14 */
  295. -      "bad starting value",            /* 15 */
  296. -      "floating exception",            /* 16 */
  297. -      "internal inconsistency (data structure)",/* 17 */
  298. -      "unexpected end-of-file",        /* 18 */
  299. -      "shared vectors (cannot release them)",  /* 19 */  
  300. -      "negative argument",            /* 20 */
  301. -      "cannot overwrite object"             /* 21 */
  302. -     };
  303. -
  304. -#define    MAXERR    (sizeof(err_mesg)/sizeof(char *))
  305. -
  306. -static char *warn_mesg[] = {
  307. -   "unknown warning",                /* 0 */
  308. -   "wrong type number (use macro TYPE_*)",    /* 1 */
  309. -   "no corresponding mem_stat_mark",        /* 2 */
  310. -   "computed norm of a residual is less than 0",  /* 3 */
  311. -   "resizing a shared vector"            /* 4 */
  312. -};
  313. -
  314. -#define MAXWARN  (sizeof(warn_mesg)/sizeof(char *))
  315. -
  316. -
  317. -
  318. -#define    MAX_ERRS    100
  319. -
  320. -jmp_buf    restart;
  321. -
  322. -
  323. -/* array of pointers to lists of errors */
  324. -
  325. -typedef struct {
  326. -   char **listp;    /* pointer to a list of errors */
  327. -   unsigned len;    /* length of the list */
  328. -   unsigned warn;   /* =FALSE - errors, =TRUE - warnings */
  329. -}  Err_list;
  330. -
  331. -static Err_list     err_list[ERR_LIST_MAX_LEN] = {
  332. - {err_mesg,MAXERR,FALSE},    /* basic errors list */
  333. - {warn_mesg,MAXWARN,TRUE}    /* basic warnings list */
  334. -};
  335. -
  336. -
  337. -static int err_list_end = 2;   /* number of elements in err_list */
  338. -
  339. -/* attach a new list of errors pointed by err_ptr
  340. -   or change a previous one;
  341. -   list_len is the number of elements in the list;
  342. -   list_num is the list number;
  343. -   warn == FALSE - errors (stop the program),
  344. -   warn == TRUE - warnings (continue the program);
  345. -   Note: lists numbered 0 and 1 are attached automatically,
  346. -   you do not need to do it
  347. -   */
  348. -int err_list_attach(list_num, list_len,err_ptr,warn)
  349. -int list_num, list_len, warn;
  350. -char **err_ptr;
  351. -{
  352. -   if (list_num < 0 || list_len <= 0 ||
  353. -       err_ptr == (char **)NULL) 
  354. -     return -1;
  355. -   
  356. -   if (list_num >= ERR_LIST_MAX_LEN) {
  357. -    fprintf(stderr,"\n file \"%s\": %s %s\n",
  358. -        "err.c","increase the value of ERR_LIST_MAX_LEN",
  359. -        "in matrix.h and zmatdef.h");
  360. -    if ( ! isatty(fileno(stdout)) )
  361. -      fprintf(stderr,"\n file \"%s\": %s %s\n",
  362. -          "err.c","increase the value of ERR_LIST_MAX_LEN",
  363. -          "in matrix.h and zmatdef.h");
  364. -    printf("Exiting program\n");
  365. -    exit(0);
  366. -     }
  367. -
  368. -   if (err_list[list_num].listp != (char **)NULL &&
  369. -       err_list[list_num].listp != err_ptr)
  370. -     free((char *)err_list[list_num].listp);
  371. -   err_list[list_num].listp = err_ptr;
  372. -   err_list[list_num].len = list_len;
  373. -   err_list[list_num].warn = warn;
  374. -   err_list_end = list_num+1;
  375. -   
  376. -   return list_num;
  377. -}
  378. -
  379. -
  380. -/* release the error list numbered list_num */
  381. -int err_list_free(list_num)
  382. -int list_num;
  383. -{
  384. -   if (list_num < 0 || list_num >= err_list_end) return -1;
  385. -   if (err_list[list_num].listp != (char **)NULL) {
  386. -      err_list[list_num].listp = (char **)NULL;
  387. -      err_list[list_num].len = 0;
  388. -      err_list[list_num].warn = 0;
  389. -   }
  390. -   return 0;
  391. -}
  392. -
  393. -
  394. -/* check if list_num is attached;
  395. -   return FALSE if not;
  396. -   return TRUE if yes
  397. -   */
  398. -int err_is_list_attached(list_num)
  399. -int list_num;
  400. -{
  401. -   if (list_num < 0 || list_num >= err_list_end)
  402. -     return FALSE;
  403. -   
  404. -   if (err_list[list_num].listp != (char **)NULL)
  405. -     return TRUE;
  406. -   
  407. -   return FALSE;
  408. -}
  409. -
  410. -/* other local variables */
  411. -
  412. -static    int    err_flag = EF_EXIT, num_errs = 0, cnt_errs = 1;
  413. -
  414. -/* set_err_flag -- sets err_flag -- returns old err_flag */
  415. -int    set_err_flag(flag)
  416. -int    flag;
  417. -{
  418. -   int    tmp;
  419. -   
  420. -   tmp = err_flag;
  421. -   err_flag = flag;
  422. -   return tmp;
  423. -}
  424. -
  425. -/* count_errs -- sets cnt_errs (TRUE/FALSE) & returns old value */
  426. -int    count_errs(flag)
  427. -int    flag;
  428. -{
  429. -   int    tmp;
  430. -   
  431. -   tmp = cnt_errs;
  432. -   cnt_errs = flag;
  433. -   return tmp;
  434. -}
  435. -
  436. -/* ev_err -- reports error (err_num) in file "file" at line "line_num" and
  437. -   returns to user error handler;
  438. -   list_num is an error list number (0 is the basic list 
  439. -   pointed by err_mesg, 1 is the basic list of warnings)
  440. - */
  441. -int    ev_err(file,err_num,line_num,fn_name,list_num)
  442. -char    *file, *fn_name;
  443. -int    err_num, line_num,list_num;
  444. -{
  445. -   int    num;
  446. -   
  447. -   if ( err_num < 0 ) err_num = 0;
  448. -   
  449. -   if (list_num < 0 || list_num >= err_list_end ||
  450. -       err_list[list_num].listp == (char **)NULL) {
  451. -      fprintf(stderr,
  452. -          "\n Not (properly) attached list of errors: list_num = %d\n",
  453. -          list_num);
  454. -      fprintf(stderr," Call \"err_list_attach\" in your program\n");
  455. -      if ( ! isatty(fileno(stdout)) ) {
  456. -     fprintf(stderr,
  457. -         "\n Not (properly) attached list of errors: list_num = %d\n",
  458. -         list_num);
  459. -     fprintf(stderr," Call \"err_list_attach\" in your program\n");
  460. -      }
  461. -      printf("\nExiting program\n");
  462. -      exit(0);
  463. -   }
  464. -   
  465. -   num = err_num;
  466. -   if ( num >= err_list[list_num].len ) num = 0;
  467. -   
  468. -   if ( cnt_errs && ++num_errs >= MAX_ERRS )    /* too many errors */
  469. -   {
  470. -      fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n",
  471. -          file,line_num,err_list[list_num].listp[num],
  472. -          isascii(*fn_name) ? fn_name : "???");
  473. -      if ( ! isatty(fileno(stdout)) )
  474. -    fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n",
  475. -        file,line_num,err_list[list_num].listp[num],
  476. -        isascii(*fn_name) ? fn_name : "???");
  477. -      printf("Sorry, too many errors: %d\n",num_errs);
  478. -      printf("Exiting program\n");
  479. -      exit(0);
  480. -   }
  481. -   if ( err_list[list_num].warn )
  482. -       switch ( err_flag )
  483. -       {
  484. -       case EF_SILENT: break;
  485. -       default:
  486. -       fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n\n",
  487. -           file,line_num,err_list[list_num].listp[num],
  488. -           isascii(*fn_name) ? fn_name : "???");
  489. -       if ( ! isatty(fileno(stdout)) )
  490. -           fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n\n",
  491. -               file,line_num,err_list[list_num].listp[num],
  492. -               isascii(*fn_name) ? fn_name : "???");
  493. -       break;
  494. -       }
  495. -   else
  496. -       switch ( err_flag )
  497. -       {
  498. -       case EF_SILENT:
  499. -       longjmp(restart,(err_num==0)? -1 : err_num);
  500. -       break;
  501. -       case EF_ABORT:
  502. -       fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n",
  503. -           file,line_num,err_list[list_num].listp[num],
  504. -           isascii(*fn_name) ? fn_name : "???");
  505. -       if ( ! isatty(fileno(stdout)) )
  506. -           fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n",
  507. -               file,line_num,err_list[list_num].listp[num],
  508. -               isascii(*fn_name) ? fn_name : "???");
  509. -       abort();
  510. -       break;
  511. -       case EF_JUMP:
  512. -       fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n",
  513. -           file,line_num,err_list[list_num].listp[num],
  514. -           isascii(*fn_name) ? fn_name : "???");
  515. -       if ( ! isatty(fileno(stdout)) )
  516. -           fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n",
  517. -               file,line_num,err_list[list_num].listp[num],
  518. -               isascii(*fn_name) ? fn_name : "???");
  519. -       longjmp(restart,(err_num==0)? -1 : err_num);
  520. -       break;
  521. -       default:
  522. -       fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n\n",
  523. -           file,line_num,err_list[list_num].listp[num],
  524. -           isascii(*fn_name) ? fn_name : "???");
  525. -       if ( ! isatty(fileno(stdout)) )
  526. -           fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n\n",
  527. -               file,line_num,err_list[list_num].listp[num],
  528. -               isascii(*fn_name) ? fn_name : "???");
  529. -       
  530. -       break;
  531. -       }
  532. -   
  533. -   /* ensure exit if fall through */
  534. -   if ( ! err_list[list_num].warn )  exit(0);
  535. -
  536. -   return 0;
  537. -}
  538. -
  539. -/* float_error -- catches floating arithmetic signals */
  540. -static void    float_error(num)
  541. -int    num;
  542. -{
  543. -   signal(SIGFPE,float_error);
  544. -   /* fprintf(stderr,"SIGFPE: signal #%d\n",num); */
  545. -   /* fprintf(stderr,"errno = %d\n",errno); */
  546. -   ev_err("???.c",E_SIGNAL,0,"???",0);
  547. -}
  548. -
  549. -/* catch_signal -- sets up float_error() to catch SIGFPE's */
  550. -void    catch_FPE()
  551. -{
  552. -   signal(SIGFPE,float_error);
  553. -}
  554. -
  555. //GO.SYSIN DD err.c
  556. echo matrixio.c 1>&2
  557. sed >matrixio.c <<'//GO.SYSIN DD matrixio.c' 's/^-//'
  558. -
  559. -/**************************************************************************
  560. -**
  561. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  562. -**
  563. -**                 Meschach Library
  564. -** 
  565. -** This Meschach Library is provided "as is" without any express 
  566. -** or implied warranty of any kind with respect to this software. 
  567. -** In particular the authors shall not be liable for any direct, 
  568. -** indirect, special, incidental or consequential damages arising 
  569. -** in any way from use of the software.
  570. -** 
  571. -** Everyone is granted permission to copy, modify and redistribute this
  572. -** Meschach Library, provided:
  573. -**  1.  All copies contain this copyright notice.
  574. -**  2.  All modified copies shall carry a notice stating who
  575. -**      made the last modification and the date of such modification.
  576. -**  3.  No charge is made for this software or works derived from it.  
  577. -**      This clause shall not be construed as constraining other software
  578. -**      distributed on the same medium as this software, nor is a
  579. -**      distribution fee considered a charge.
  580. -**
  581. -***************************************************************************/
  582. -
  583. -
  584. -/* 1.6 matrixio.c 11/25/87 */
  585. -
  586. -
  587. -#include        <stdio.h>
  588. -#include        <ctype.h>
  589. -#include        "matrix.h"
  590. -
  591. -static char rcsid[] = "$Id: matrixio.c,v 1.4 1994/01/13 05:31:10 des Exp $";
  592. -
  593. -
  594. -/* local variables */
  595. -static char line[MAXLINE];
  596. -
  597. -
  598. -/**************************************************************************
  599. -  Input routines
  600. -  **************************************************************************/
  601. -/* skipjunk -- skips white spaces and strings of the form #....\n
  602. -   Here .... is a comment string */
  603. -int     skipjunk(fp)
  604. -FILE    *fp;
  605. -{
  606. -     int        c;
  607. -     
  608. -     for ( ; ; )        /* forever do... */
  609. -     {
  610. -      /* skip blanks */
  611. -      do
  612. -           c = getc(fp);
  613. -      while ( isspace(c) );
  614. -      
  615. -      /* skip comments (if any) */
  616. -      if ( c == '#' )
  617. -           /* yes it is a comment (line) */
  618. -           while ( (c=getc(fp)) != '\n' )
  619. -            ;
  620. -      else
  621. -      {
  622. -           ungetc(c,fp);
  623. -           break;
  624. -      }
  625. -     }
  626. -     return 0;
  627. -}
  628. -
  629. -MAT     *m_finput(fp,a)
  630. -FILE    *fp;
  631. -MAT     *a;
  632. -{
  633. -     MAT        *im_finput(),*bm_finput();
  634. -     
  635. -     if ( isatty(fileno(fp)) )
  636. -      return im_finput(fp,a);
  637. -     else
  638. -      return bm_finput(fp,a);
  639. -}
  640. -
  641. -/* im_finput -- interactive input of matrix */
  642. -MAT     *im_finput(fp,mat)
  643. -FILE    *fp;
  644. -MAT     *mat;
  645. -{
  646. -     char       c;
  647. -     u_int      i, j, m, n, dynamic;
  648. -     /* dynamic set to TRUE if memory allocated here */
  649. -     
  650. -     /* get matrix size */
  651. -     if ( mat != (MAT *)NULL && mat->m<MAXDIM && mat->n<MAXDIM )
  652. -     {  m = mat->m;     n = mat->n;     dynamic = FALSE;        }
  653. -     else
  654. -     {
  655. -      dynamic = TRUE;
  656. -      do
  657. -      {
  658. -           fprintf(stderr,"Matrix: rows cols:");
  659. -           if ( fgets(line,MAXLINE,fp)==NULL )
  660. -            error(E_INPUT,"im_finput");
  661. -      } while ( sscanf(line,"%u%u",&m,&n)<2 || m>MAXDIM || n>MAXDIM );
  662. -      mat = m_get(m,n);
  663. -     }
  664. -     
  665. -     /* input elements */
  666. -     for ( i=0; i<m; i++ )
  667. -     {
  668. -     redo:
  669. -      fprintf(stderr,"row %u:\n",i);
  670. -      for ( j=0; j<n; j++ )
  671. -           do
  672. -           {
  673. -           redo2:
  674. -            fprintf(stderr,"entry (%u,%u): ",i,j);
  675. -            if ( !dynamic )
  676. -             fprintf(stderr,"old %14.9g new: ",
  677. -                 mat->me[i][j]);
  678. -            if ( fgets(line,MAXLINE,fp)==NULL )
  679. -             error(E_INPUT,"im_finput");
  680. -            if ( (*line == 'b' || *line == 'B') && j > 0 )
  681. -            {   j--;    dynamic = FALSE;        goto redo2;     }
  682. -            if ( (*line == 'f' || *line == 'F') && j < n-1 )
  683. -            {   j++;    dynamic = FALSE;        goto redo2;     }
  684. -#if REAL == DOUBLE
  685. -           } while ( *line=='\0' || sscanf(line,"%lf",&mat->me[i][j])<1 );
  686. -#elif REAL == FLOAT
  687. -           } while ( *line=='\0' || sscanf(line,"%f",&mat->me[i][j])<1 );
  688. -#endif
  689. -      fprintf(stderr,"Continue: ");
  690. -      fscanf(fp,"%c",&c);
  691. -      if ( c == 'n' || c == 'N' )
  692. -      {    dynamic = FALSE;                 goto redo;      }
  693. -      if ( (c == 'b' || c == 'B') /* && i > 0 */ )
  694. -      {     if ( i > 0 )
  695. -            i--;
  696. -        dynamic = FALSE;        goto redo;
  697. -      }
  698. -     }
  699. -     
  700. -     return (mat);
  701. -}
  702. -
  703. -/* bm_finput -- batch-file input of matrix */
  704. -MAT     *bm_finput(fp,mat)
  705. -FILE    *fp;
  706. -MAT     *mat;
  707. -{
  708. -     u_int      i,j,m,n,dummy;
  709. -     int        io_code;
  710. -     
  711. -     /* get dimension */
  712. -     skipjunk(fp);
  713. -     if ((io_code=fscanf(fp," Matrix: %u by %u",&m,&n)) < 2 ||
  714. -     m>MAXDIM || n>MAXDIM )
  715. -      error(io_code==EOF ? E_EOF : E_FORMAT,"bm_finput");
  716. -     
  717. -     /* allocate memory if necessary */
  718. -     if ( mat==(MAT *)NULL )
  719. -      mat = m_resize(mat,m,n);
  720. -     
  721. -     /* get entries */
  722. -     for ( i=0; i<m; i++ )
  723. -     {
  724. -      skipjunk(fp);
  725. -      if ( fscanf(fp," row %u:",&dummy) < 1 )
  726. -           error(E_FORMAT,"bm_finput");
  727. -      for ( j=0; j<n; j++ )
  728. -#if REAL == DOUBLE
  729. -           if ((io_code=fscanf(fp,"%lf",&mat->me[i][j])) < 1 )
  730. -#elif REAL == FLOAT
  731. -           if ((io_code=fscanf(fp,"%f",&mat->me[i][j])) < 1 )
  732. -#endif
  733. -            error(io_code==EOF ? 7 : 6,"bm_finput");
  734. -     }
  735. -     
  736. -     return (mat);
  737. -}
  738. -
  739. -PERM    *px_finput(fp,px)
  740. -FILE    *fp;
  741. -PERM    *px;
  742. -{
  743. -     PERM       *ipx_finput(),*bpx_finput();
  744. -     
  745. -     if ( isatty(fileno(fp)) )
  746. -      return ipx_finput(fp,px);
  747. -     else
  748. -      return bpx_finput(fp,px);
  749. -}
  750. -
  751. -
  752. -/* ipx_finput -- interactive input of permutation */
  753. -PERM    *ipx_finput(fp,px)
  754. -FILE    *fp;
  755. -PERM    *px;
  756. -{
  757. -     u_int      i,j,size,dynamic; /* dynamic set if memory allocated here */
  758. -     u_int      entry,ok;
  759. -     
  760. -     /* get permutation size */
  761. -     if ( px!=(PERM *)NULL && px->size<MAXDIM )
  762. -     {  size = px->size;        dynamic = FALSE;        }
  763. -     else
  764. -     {
  765. -      dynamic = TRUE;
  766. -      do
  767. -      {
  768. -           fprintf(stderr,"Permutation: size: ");
  769. -           if ( fgets(line,MAXLINE,fp)==NULL )
  770. -            error(E_INPUT,"ipx_finput");
  771. -      } while ( sscanf(line,"%u",&size)<1 || size>MAXDIM );
  772. -      px = px_get(size);
  773. -     }
  774. -     
  775. -     /* get entries */
  776. -     i = 0;
  777. -     while ( i<size )
  778. -     {
  779. -      /* input entry */
  780. -      do
  781. -      {
  782. -      redo:
  783. -           fprintf(stderr,"entry %u: ",i);
  784. -           if ( !dynamic )
  785. -            fprintf(stderr,"old: %u->%u new: ",
  786. -                i,px->pe[i]);
  787. -           if ( fgets(line,MAXLINE,fp)==NULL )
  788. -            error(E_INPUT,"ipx_finput");
  789. -           if ( (*line == 'b' || *line == 'B') && i > 0 )
  790. -           {        i--;    dynamic = FALSE;        goto redo;      }
  791. -      } while ( *line=='\0' || sscanf(line,"%u",&entry) < 1 );
  792. -      /* check entry */
  793. -      ok = (entry < size);
  794. -      for ( j=0; j<i; j++ )
  795. -           ok &= (entry != px->pe[j]);
  796. -      if ( ok )
  797. -      {
  798. -           px->pe[i] = entry;
  799. -           i++;
  800. -      }
  801. -     }
  802. -     
  803. -     return (px);
  804. -}
  805. -
  806. -/* bpx_finput -- batch-file input of permutation */
  807. -PERM    *bpx_finput(fp,px)
  808. -FILE    *fp;
  809. -PERM    *px;
  810. -{
  811. -     u_int      i,j,size,entry,ok;
  812. -     int        io_code;
  813. -     
  814. -     /* get size of permutation */
  815. -     skipjunk(fp);
  816. -     if ((io_code=fscanf(fp," Permutation: size:%u",&size)) < 1 ||
  817. -     size>MAXDIM )
  818. -      error(io_code==EOF ? 7 : 6,"bpx_finput");
  819. -     
  820. -     /* allocate memory if necessary */
  821. -     if ( px==(PERM *)NULL || px->size<size )
  822. -      px = px_resize(px,size);
  823. -     
  824. -     /* get entries */
  825. -     skipjunk(fp);
  826. -     i = 0;
  827. -     while ( i<size )
  828. -     {
  829. -      /* input entry */
  830. -      if ((io_code=fscanf(fp,"%*u -> %u",&entry)) < 1 )
  831. -           error(io_code==EOF ? 7 : 6,"bpx_finput");
  832. -      /* check entry */
  833. -      ok = (entry < size);
  834. -      for ( j=0; j<i; j++ )
  835. -           ok &= (entry != px->pe[j]);
  836. -      if ( ok )
  837. -      {
  838. -           px->pe[i] = entry;
  839. -           i++;
  840. -      }
  841. -      else
  842. -           error(E_BOUNDS,"bpx_finput");
  843. -     }
  844. -     
  845. -     return (px);
  846. -}
  847. -
  848. -
  849. -VEC     *v_finput(fp,x)
  850. -FILE    *fp;
  851. -VEC     *x;
  852. -{
  853. -     VEC        *ifin_vec(),*bfin_vec();
  854. -     
  855. -     if ( isatty(fileno(fp)) )
  856. -      return ifin_vec(fp,x);
  857. -     else
  858. -      return bfin_vec(fp,x);
  859. -}
  860. -
  861. -/* ifin_vec -- interactive input of vector */
  862. -VEC     *ifin_vec(fp,vec)
  863. -FILE    *fp;
  864. -VEC     *vec;
  865. -{
  866. -     u_int      i,dim,dynamic;  /* dynamic set if memory allocated here */
  867. -     
  868. -     /* get vector dimension */
  869. -     if ( vec != (VEC *)NULL && vec->dim<MAXDIM )
  870. -     {  dim = vec->dim; dynamic = FALSE;        }
  871. -     else
  872. -     {
  873. -      dynamic = TRUE;
  874. -      do
  875. -      {
  876. -           fprintf(stderr,"Vector: dim: ");
  877. -           if ( fgets(line,MAXLINE,fp)==NULL )
  878. -            error(E_INPUT,"ifin_vec");
  879. -      } while ( sscanf(line,"%u",&dim)<1 || dim>MAXDIM );
  880. -      vec = v_get(dim);
  881. -     }
  882. -     
  883. -     /* input elements */
  884. -     for ( i=0; i<dim; i++ )
  885. -      do
  886. -      {
  887. -      redo:
  888. -           fprintf(stderr,"entry %u: ",i);
  889. -           if ( !dynamic )
  890. -            fprintf(stderr,"old %14.9g new: ",vec->ve[i]);
  891. -           if ( fgets(line,MAXLINE,fp)==NULL )
  892. -            error(E_INPUT,"ifin_vec");
  893. -           if ( (*line == 'b' || *line == 'B') && i > 0 )
  894. -           {        i--;    dynamic = FALSE;        goto redo;         }
  895. -           if ( (*line == 'f' || *line == 'F') && i < dim-1 )
  896. -           {        i++;    dynamic = FALSE;        goto redo;         }
  897. -#if REAL == DOUBLE
  898. -      } while ( *line=='\0' || sscanf(line,"%lf",&vec->ve[i]) < 1 );
  899. -#elif REAL == FLOAT
  900. -          } while ( *line=='\0' || sscanf(line,"%f",&vec->ve[i]) < 1 );
  901. -#endif
  902. -     
  903. -     return (vec);
  904. -}
  905. -
  906. -/* bfin_vec -- batch-file input of vector */
  907. -VEC     *bfin_vec(fp,vec)
  908. -FILE    *fp;
  909. -VEC     *vec;
  910. -{
  911. -     u_int      i,dim;
  912. -     int        io_code;
  913. -     
  914. -     /* get dimension */
  915. -     skipjunk(fp);
  916. -     if ((io_code=fscanf(fp," Vector: dim:%u",&dim)) < 1 ||
  917. -     dim>MAXDIM )
  918. -      error(io_code==EOF ? 7 : 6,"bfin_vec");
  919. -     
  920. -     /* allocate memory if necessary */
  921. -     if ( vec==(VEC *)NULL )
  922. -      vec = v_resize(vec,dim);
  923. -     
  924. -     /* get entries */
  925. -     skipjunk(fp);
  926. -     for ( i=0; i<dim; i++ )
  927. -#if REAL == DOUBLE
  928. -      if ((io_code=fscanf(fp,"%lf",&vec->ve[i])) < 1 )
  929. -#elif REAL == FLOAT
  930. -      if ((io_code=fscanf(fp,"%f",&vec->ve[i])) < 1 )
  931. -#endif
  932. -           error(io_code==EOF ? 7 : 6,"bfin_vec");
  933. -     
  934. -     return (vec);
  935. -}
  936. -
  937. -/**************************************************************************
  938. -  Output routines
  939. -  **************************************************************************/
  940. -static char    *format = "%14.9g ";
  941. -
  942. -char    *setformat(f_string)
  943. -char    *f_string;
  944. -{
  945. -    char    *old_f_string;
  946. -    old_f_string = format;
  947. -    if ( f_string != (char *)NULL && *f_string != '\0' )
  948. -    format = f_string;
  949. -
  950. -    return old_f_string;
  951. -}
  952. -
  953. -void    m_foutput(fp,a)
  954. -FILE    *fp;
  955. -MAT     *a;
  956. -{
  957. -     u_int      i, j, tmp;
  958. -     
  959. -     if ( a == (MAT *)NULL )
  960. -     {  fprintf(fp,"Matrix: NULL\n");   return;         }
  961. -     fprintf(fp,"Matrix: %d by %d\n",a->m,a->n);
  962. -     if ( a->me == (Real **)NULL )
  963. -     {  fprintf(fp,"NULL\n");           return;         }
  964. -     for ( i=0; i<a->m; i++ )   /* for each row... */
  965. -     {
  966. -      fprintf(fp,"row %u: ",i);
  967. -      for ( j=0, tmp=2; j<a->n; j++, tmp++ )
  968. -      {             /* for each col in row... */
  969. -           fprintf(fp,format,a->me[i][j]);
  970. -           if ( ! (tmp % 5) )       putc('\n',fp);
  971. -      }
  972. -      if ( tmp % 5 != 1 )   putc('\n',fp);
  973. -     }
  974. -}
  975. -
  976. -void    px_foutput(fp,px)
  977. -FILE    *fp;
  978. -PERM    *px;
  979. -{
  980. -     u_int      i;
  981. -     
  982. -     if ( px == (PERM *)NULL )
  983. -     {  fprintf(fp,"Permutation: NULL\n");      return;         }
  984. -     fprintf(fp,"Permutation: size: %u\n",px->size);
  985. -     if ( px->pe == (u_int *)NULL )
  986. -     {  fprintf(fp,"NULL\n");   return;         }
  987. -     for ( i=0; i<px->size; i++ )
  988. -    if ( ! (i % 8) && i != 0 )
  989. -      fprintf(fp,"\n  %u->%u ",i,px->pe[i]);
  990. -    else
  991. -      fprintf(fp,"%u->%u ",i,px->pe[i]);
  992. -     fprintf(fp,"\n");
  993. -}
  994. -
  995. -void    v_foutput(fp,x)
  996. -FILE    *fp;
  997. -VEC     *x;
  998. -{
  999. -     u_int      i, tmp;
  1000. -     
  1001. -     if ( x == (VEC *)NULL )
  1002. -     {  fprintf(fp,"Vector: NULL\n");   return;         }
  1003. -     fprintf(fp,"Vector: dim: %d\n",x->dim);
  1004. -     if ( x->ve == (Real *)NULL )
  1005. -     {  fprintf(fp,"NULL\n");   return;         }
  1006. -     for ( i=0, tmp=0; i<x->dim; i++, tmp++ )
  1007. -     {
  1008. -      fprintf(fp,format,x->ve[i]);
  1009. -      if ( tmp % 5 == 4 )   putc('\n',fp);
  1010. -     }
  1011. -     if ( tmp % 5 != 0 )        putc('\n',fp);
  1012. -}
  1013. -
  1014. -
  1015. -void    m_dump(fp,a)
  1016. -FILE    *fp;
  1017. -MAT     *a;
  1018. -{
  1019. -    u_int   i, j, tmp;
  1020. -     
  1021. -     if ( a == (MAT *)NULL )
  1022. -     {  fprintf(fp,"Matrix: NULL\n");   return;         }
  1023. -     fprintf(fp,"Matrix: %d by %d @ 0x%lx\n",a->m,a->n,(long)a);
  1024. -     fprintf(fp,"\tmax_m = %d, max_n = %d, max_size = %d\n",
  1025. -         a->max_m, a->max_n, a->max_size);
  1026. -     if ( a->me == (Real **)NULL )
  1027. -     {  fprintf(fp,"NULL\n");           return;         }
  1028. -     fprintf(fp,"a->me @ 0x%lx\n",(long)(a->me));
  1029. -     fprintf(fp,"a->base @ 0x%lx\n",(long)(a->base));
  1030. -     for ( i=0; i<a->m; i++ )   /* for each row... */
  1031. -     {
  1032. -      fprintf(fp,"row %u: @ 0x%lx ",i,(long)(a->me[i]));
  1033. -      for ( j=0, tmp=2; j<a->n; j++, tmp++ )
  1034. -      {             /* for each col in row... */
  1035. -           fprintf(fp,format,a->me[i][j]);
  1036. -           if ( ! (tmp % 5) )       putc('\n',fp);
  1037. -      }
  1038. -      if ( tmp % 5 != 1 )   putc('\n',fp);
  1039. -     }
  1040. -}
  1041. -
  1042. -void    px_dump(fp,px)
  1043. -FILE    *fp;
  1044. -PERM    *px;
  1045. -{
  1046. -     u_int      i;
  1047. -     
  1048. -     if ( ! px )
  1049. -     {  fprintf(fp,"Permutation: NULL\n");      return;         }
  1050. -     fprintf(fp,"Permutation: size: %u @ 0x%lx\n",px->size,(long)(px));
  1051. -     if ( ! px->pe )
  1052. -     {  fprintf(fp,"NULL\n");   return;         }
  1053. -     fprintf(fp,"px->pe @ 0x%lx\n",(long)(px->pe));
  1054. -     for ( i=0; i<px->size; i++ )
  1055. -      fprintf(fp,"%u->%u ",i,px->pe[i]);
  1056. -     fprintf(fp,"\n");
  1057. -}
  1058. -
  1059. -
  1060. -void    v_dump(fp,x)
  1061. -FILE    *fp;
  1062. -VEC     *x;
  1063. -{
  1064. -     u_int      i, tmp;
  1065. -     
  1066. -     if ( ! x )
  1067. -     {  fprintf(fp,"Vector: NULL\n");   return;         }
  1068. -     fprintf(fp,"Vector: dim: %d @ 0x%lx\n",x->dim,(long)(x));
  1069. -     if ( ! x->ve )
  1070. -     {  fprintf(fp,"NULL\n");   return;         }
  1071. -     fprintf(fp,"x->ve @ 0x%lx\n",(long)(x->ve));
  1072. -     for ( i=0, tmp=0; i<x->dim; i++, tmp++ )
  1073. -     {
  1074. -      fprintf(fp,format,x->ve[i]);
  1075. -      if ( tmp % 5 == 4 )   putc('\n',fp);
  1076. -     }
  1077. -     if ( tmp % 5 != 0 )        putc('\n',fp);
  1078. -}
  1079. -
  1080. //GO.SYSIN DD matrixio.c
  1081. echo memory.c 1>&2
  1082. sed >memory.c <<'//GO.SYSIN DD memory.c' 's/^-//'
  1083. -
  1084. -/**************************************************************************
  1085. -**
  1086. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1087. -**
  1088. -**                 Meschach Library
  1089. -** 
  1090. -** This Meschach Library is provided "as is" without any express 
  1091. -** or implied warranty of any kind with respect to this software. 
  1092. -** In particular the authors shall not be liable for any direct, 
  1093. -** indirect, special, incidental or consequential damages arising 
  1094. -** in any way from use of the software.
  1095. -** 
  1096. -** Everyone is granted permission to copy, modify and redistribute this
  1097. -** Meschach Library, provided:
  1098. -**  1.  All copies contain this copyright notice.
  1099. -**  2.  All modified copies shall carry a notice stating who
  1100. -**      made the last modification and the date of such modification.
  1101. -**  3.  No charge is made for this software or works derived from it.  
  1102. -**      This clause shall not be construed as constraining other software
  1103. -**      distributed on the same medium as this software, nor is a
  1104. -**      distribution fee considered a charge.
  1105. -**
  1106. -***************************************************************************/
  1107. -
  1108. -
  1109. -/* memory.c 1.3 11/25/87 */
  1110. -
  1111. -#include     "matrix.h"
  1112. -
  1113. -
  1114. -static    char    rcsid[] = "$Id: memory.c,v 1.13 1994/04/05 02:10:37 des Exp $";
  1115. -
  1116. -/* m_get -- gets an mxn matrix (in MAT form) by dynamic memory allocation */
  1117. -MAT    *m_get(m,n)
  1118. -int    m,n;
  1119. -{
  1120. -   MAT    *matrix;
  1121. -   int    i;
  1122. -   
  1123. -   if (m < 0 || n < 0)
  1124. -     error(E_NEG,"m_get");
  1125. -
  1126. -   if ((matrix=NEW(MAT)) == (MAT *)NULL )
  1127. -     error(E_MEM,"m_get");
  1128. -   else if (mem_info_is_on()) {
  1129. -      mem_bytes(TYPE_MAT,0,sizeof(MAT));
  1130. -      mem_numvar(TYPE_MAT,1);
  1131. -   }
  1132. -   
  1133. -   matrix->m = m;        matrix->n = matrix->max_n = n;
  1134. -   matrix->max_m = m;    matrix->max_size = m*n;
  1135. -#ifndef SEGMENTED
  1136. -   if ((matrix->base = NEW_A(m*n,Real)) == (Real *)NULL )
  1137. -   {
  1138. -      free(matrix);
  1139. -      error(E_MEM,"m_get");
  1140. -   }
  1141. -   else if (mem_info_is_on()) {
  1142. -      mem_bytes(TYPE_MAT,0,m*n*sizeof(Real));
  1143. -   }
  1144. -#else
  1145. -   matrix->base = (Real *)NULL;
  1146. -#endif
  1147. -   if ((matrix->me = (Real **)calloc(m,sizeof(Real *))) == 
  1148. -       (Real **)NULL )
  1149. -   {    free(matrix->base);    free(matrix);
  1150. -    error(E_MEM,"m_get");
  1151. -     }
  1152. -   else if (mem_info_is_on()) {
  1153. -      mem_bytes(TYPE_MAT,0,m*sizeof(Real *));
  1154. -   }
  1155. -   
  1156. -#ifndef SEGMENTED
  1157. -   /* set up pointers */
  1158. -   for ( i=0; i<m; i++ )
  1159. -     matrix->me[i] = &(matrix->base[i*n]);
  1160. -#else
  1161. -   for ( i = 0; i < m; i++ )
  1162. -     if ( (matrix->me[i]=NEW_A(n,Real)) == (Real *)NULL )
  1163. -       error(E_MEM,"m_get");
  1164. -     else if (mem_info_is_on()) {
  1165. -    mem_bytes(TYPE_MAT,0,n*sizeof(Real));
  1166. -       }
  1167. -#endif
  1168. -   
  1169. -   return (matrix);
  1170. -}
  1171. -
  1172. -
  1173. -/* px_get -- gets a PERM of given 'size' by dynamic memory allocation
  1174. -   -- Note: initialized to the identity permutation */
  1175. -PERM    *px_get(size)
  1176. -int    size;
  1177. -{
  1178. -   PERM    *permute;
  1179. -   int    i;
  1180. -
  1181. -   if (size < 0)
  1182. -     error(E_NEG,"px_get");
  1183. -
  1184. -   if ((permute=NEW(PERM)) == (PERM *)NULL )
  1185. -     error(E_MEM,"px_get");
  1186. -   else if (mem_info_is_on()) {
  1187. -      mem_bytes(TYPE_PERM,0,sizeof(PERM));
  1188. -      mem_numvar(TYPE_PERM,1);
  1189. -   }
  1190. -   
  1191. -   permute->size = permute->max_size = size;
  1192. -   if ((permute->pe = NEW_A(size,u_int)) == (u_int *)NULL )
  1193. -     error(E_MEM,"px_get");
  1194. -   else if (mem_info_is_on()) {
  1195. -      mem_bytes(TYPE_PERM,0,size*sizeof(u_int));
  1196. -   }
  1197. -   
  1198. -   for ( i=0; i<size; i++ )
  1199. -     permute->pe[i] = i;
  1200. -   
  1201. -   return (permute);
  1202. -}
  1203. -
  1204. -/* v_get -- gets a VEC of dimension 'dim'
  1205. -   -- Note: initialized to zero */
  1206. -VEC    *v_get(size)
  1207. -int    size;
  1208. -{
  1209. -   VEC    *vector;
  1210. -   
  1211. -   if (size < 0)
  1212. -     error(E_NEG,"v_get");
  1213. -
  1214. -   if ((vector=NEW(VEC)) == (VEC *)NULL )
  1215. -     error(E_MEM,"v_get");
  1216. -   else if (mem_info_is_on()) {
  1217. -      mem_bytes(TYPE_VEC,0,sizeof(VEC));
  1218. -      mem_numvar(TYPE_VEC,1);
  1219. -   }
  1220. -   
  1221. -   vector->dim = vector->max_dim = size;
  1222. -   if ((vector->ve=NEW_A(size,Real)) == (Real *)NULL )
  1223. -   {
  1224. -      free(vector);
  1225. -      error(E_MEM,"v_get");
  1226. -   }
  1227. -   else if (mem_info_is_on()) {
  1228. -      mem_bytes(TYPE_VEC,0,size*sizeof(Real));
  1229. -   }
  1230. -   
  1231. -   return (vector);
  1232. -}
  1233. -
  1234. -/* m_free -- returns MAT & asoociated memory back to memory heap */
  1235. -int    m_free(mat)
  1236. -MAT    *mat;
  1237. -{
  1238. -#ifdef SEGMENTED
  1239. -   int    i;
  1240. -#endif
  1241. -   
  1242. -   if ( mat==(MAT *)NULL || (int)(mat->m) < 0 ||
  1243. -       (int)(mat->n) < 0 )
  1244. -     /* don't trust it */
  1245. -     return (-1);
  1246. -   
  1247. -#ifndef SEGMENTED
  1248. -   if ( mat->base != (Real *)NULL ) {
  1249. -      if (mem_info_is_on()) {
  1250. -     mem_bytes(TYPE_MAT,mat->max_m*mat->max_n*sizeof(Real),0);
  1251. -      }
  1252. -      free((char *)(mat->base));
  1253. -   }
  1254. -#else
  1255. -   for ( i = 0; i < mat->max_m; i++ )
  1256. -     if ( mat->me[i] != (Real *)NULL ) {
  1257. -    if (mem_info_is_on()) {
  1258. -       mem_bytes(TYPE_MAT,mat->max_n*sizeof(Real),0);
  1259. -    }
  1260. -    free((char *)(mat->me[i]));
  1261. -     }
  1262. -#endif
  1263. -   if ( mat->me != (Real **)NULL ) {
  1264. -      if (mem_info_is_on()) {
  1265. -     mem_bytes(TYPE_MAT,mat->max_m*sizeof(Real *),0);
  1266. -      }
  1267. -      free((char *)(mat->me));
  1268. -   }
  1269. -   
  1270. -   if (mem_info_is_on()) {
  1271. -      mem_bytes(TYPE_MAT,sizeof(MAT),0);
  1272. -      mem_numvar(TYPE_MAT,-1);
  1273. -   }
  1274. -   free((char *)mat);
  1275. -   
  1276. -   return (0);
  1277. -}
  1278. -
  1279. -
  1280. -
  1281. -/* px_free -- returns PERM & asoociated memory back to memory heap */
  1282. -int    px_free(px)
  1283. -PERM    *px;
  1284. -{
  1285. -   if ( px==(PERM *)NULL || (int)(px->size) < 0 )
  1286. -     /* don't trust it */
  1287. -     return (-1);
  1288. -   
  1289. -   if ( px->pe == (u_int *)NULL ) {
  1290. -      if (mem_info_is_on()) {
  1291. -     mem_bytes(TYPE_PERM,sizeof(PERM),0);
  1292. -     mem_numvar(TYPE_PERM,-1);
  1293. -      }      
  1294. -      free((char *)px);
  1295. -   }
  1296. -   else
  1297. -   {
  1298. -      if (mem_info_is_on()) {
  1299. -     mem_bytes(TYPE_PERM,sizeof(PERM)+px->max_size*sizeof(u_int),0);
  1300. -     mem_numvar(TYPE_PERM,-1);
  1301. -      }
  1302. -      free((char *)px->pe);
  1303. -      free((char *)px);
  1304. -   }
  1305. -   
  1306. -   return (0);
  1307. -}
  1308. -
  1309. -
  1310. -
  1311. -/* v_free -- returns VEC & asoociated memory back to memory heap */
  1312. -int    v_free(vec)
  1313. -VEC    *vec;
  1314. -{
  1315. -   if ( vec==(VEC *)NULL || (int)(vec->dim) < 0 )
  1316. -     /* don't trust it */
  1317. -     return (-1);
  1318. -   
  1319. -   if ( vec->ve == (Real *)NULL ) {
  1320. -      if (mem_info_is_on()) {
  1321. -     mem_bytes(TYPE_VEC,sizeof(VEC),0);
  1322. -     mem_numvar(TYPE_VEC,-1);
  1323. -      }
  1324. -      free((char *)vec);
  1325. -   }
  1326. -   else
  1327. -   {
  1328. -      if (mem_info_is_on()) {
  1329. -     mem_bytes(TYPE_VEC,sizeof(VEC)+vec->max_dim*sizeof(Real),0);
  1330. -     mem_numvar(TYPE_VEC,-1);
  1331. -      }
  1332. -      free((char *)vec->ve);
  1333. -      free((char *)vec);
  1334. -   }
  1335. -   
  1336. -   return (0);
  1337. -}
  1338. -
  1339. -
  1340. -
  1341. -/* m_resize -- returns the matrix A of size new_m x new_n; A is zeroed
  1342. -   -- if A == NULL on entry then the effect is equivalent to m_get() */
  1343. -MAT    *m_resize(A,new_m,new_n)
  1344. -MAT    *A;
  1345. -int    new_m, new_n;
  1346. -{
  1347. -   int    i;
  1348. -   int    new_max_m, new_max_n, new_size, old_m, old_n;
  1349. -   
  1350. -   if (new_m < 0 || new_n < 0)
  1351. -     error(E_NEG,"m_resize");
  1352. -
  1353. -   if ( ! A )
  1354. -     return m_get(new_m,new_n);
  1355. -
  1356. -   /* nothing was changed */
  1357. -   if (new_m == A->m && new_n == A->n)
  1358. -     return A;
  1359. -
  1360. -   old_m = A->m;    old_n = A->n;
  1361. -   if ( new_m > A->max_m )
  1362. -   {    /* re-allocate A->me */
  1363. -      if (mem_info_is_on()) {
  1364. -     mem_bytes(TYPE_MAT,A->max_m*sizeof(Real *),
  1365. -              new_m*sizeof(Real *));
  1366. -      }
  1367. -
  1368. -      A->me = RENEW(A->me,new_m,Real *);
  1369. -      if ( ! A->me )
  1370. -    error(E_MEM,"m_resize");
  1371. -   }
  1372. -   new_max_m = max(new_m,A->max_m);
  1373. -   new_max_n = max(new_n,A->max_n);
  1374. -   
  1375. -#ifndef SEGMENTED
  1376. -   new_size = new_max_m*new_max_n;
  1377. -   if ( new_size > A->max_size )
  1378. -   {    /* re-allocate A->base */
  1379. -      if (mem_info_is_on()) {
  1380. -     mem_bytes(TYPE_MAT,A->max_m*A->max_n*sizeof(Real),
  1381. -              new_size*sizeof(Real));
  1382. -      }
  1383. -
  1384. -      A->base = RENEW(A->base,new_size,Real);
  1385. -      if ( ! A->base )
  1386. -    error(E_MEM,"m_resize");
  1387. -      A->max_size = new_size;
  1388. -   }
  1389. -   
  1390. -   /* now set up A->me[i] */
  1391. -   for ( i = 0; i < new_m; i++ )
  1392. -     A->me[i] = &(A->base[i*new_n]);
  1393. -   
  1394. -   /* now shift data in matrix */
  1395. -   if ( old_n > new_n )
  1396. -   {
  1397. -      for ( i = 1; i < min(old_m,new_m); i++ )
  1398. -    MEM_COPY((char *)&(A->base[i*old_n]),
  1399. -         (char *)&(A->base[i*new_n]),
  1400. -         sizeof(Real)*new_n);
  1401. -   }
  1402. -   else if ( old_n < new_n )
  1403. -   {
  1404. -      for ( i = (int)(min(old_m,new_m))-1; i > 0; i-- )
  1405. -      {   /* copy & then zero extra space */
  1406. -     MEM_COPY((char *)&(A->base[i*old_n]),
  1407. -          (char *)&(A->base[i*new_n]),
  1408. -          sizeof(Real)*old_n);
  1409. -     __zero__(&(A->base[i*new_n+old_n]),(new_n-old_n));
  1410. -      }
  1411. -      __zero__(&(A->base[old_n]),(new_n-old_n));
  1412. -      A->max_n = new_n;
  1413. -   }
  1414. -   /* zero out the new rows.. */
  1415. -   for ( i = old_m; i < new_m; i++ )
  1416. -     __zero__(&(A->base[i*new_n]),new_n);
  1417. -#else
  1418. -   if ( A->max_n < new_n )
  1419. -   {
  1420. -      Real    *tmp;
  1421. -      
  1422. -      for ( i = 0; i < A->max_m; i++ )
  1423. -      {
  1424. -     if (mem_info_is_on()) {
  1425. -        mem_bytes(TYPE_MAT,A->max_n*sizeof(Real),
  1426. -             new_max_n*sizeof(Real));
  1427. -     }    
  1428. -
  1429. -     if ( (tmp = RENEW(A->me[i],new_max_n,Real)) == NULL )
  1430. -       error(E_MEM,"m_resize");
  1431. -     else {    
  1432. -        A->me[i] = tmp;
  1433. -     }
  1434. -      }
  1435. -      for ( i = A->max_m; i < new_max_m; i++ )
  1436. -      {
  1437. -     if ( (tmp = NEW_A(new_max_n,Real)) == NULL )
  1438. -       error(E_MEM,"m_resize");
  1439. -     else {
  1440. -        A->me[i] = tmp;
  1441. -
  1442. -        if (mem_info_is_on()) {
  1443. -           mem_bytes(TYPE_MAT,0,new_max_n*sizeof(Real));
  1444. -        }        
  1445. -     }
  1446. -      }
  1447. -   }
  1448. -   else if ( A->max_m < new_m )
  1449. -   {
  1450. -      for ( i = A->max_m; i < new_m; i++ ) 
  1451. -    if ( (A->me[i] = NEW_A(new_max_n,Real)) == NULL )
  1452. -      error(E_MEM,"m_resize");
  1453. -    else if (mem_info_is_on()) {
  1454. -       mem_bytes(TYPE_MAT,0,new_max_n*sizeof(Real));
  1455. -    }
  1456. -      
  1457. -   }
  1458. -   
  1459. -   if ( old_n < new_n )
  1460. -   {
  1461. -      for ( i = 0; i < old_m; i++ )
  1462. -    __zero__(&(A->me[i][old_n]),new_n-old_n);
  1463. -   }
  1464. -   
  1465. -   /* zero out the new rows.. */
  1466. -   for ( i = old_m; i < new_m; i++ )
  1467. -     __zero__(A->me[i],new_n);
  1468. -#endif
  1469. -   
  1470. -   A->max_m = new_max_m;
  1471. -   A->max_n = new_max_n;
  1472. -   A->max_size = A->max_m*A->max_n;
  1473. -   A->m = new_m;    A->n = new_n;
  1474. -   
  1475. -   return A;
  1476. -}
  1477. -
  1478. -/* px_resize -- returns the permutation px with size new_size
  1479. -   -- px is set to the identity permutation */
  1480. -PERM    *px_resize(px,new_size)
  1481. -PERM    *px;
  1482. -int    new_size;
  1483. -{
  1484. -   int    i;
  1485. -   
  1486. -   if (new_size < 0)
  1487. -     error(E_NEG,"px_resize");
  1488. -
  1489. -   if ( ! px )
  1490. -     return px_get(new_size);
  1491. -   
  1492. -   /* nothing is changed */
  1493. -   if (new_size == px->size)
  1494. -     return px;
  1495. -
  1496. -   if ( new_size > px->max_size )
  1497. -   {
  1498. -      if (mem_info_is_on()) {
  1499. -     mem_bytes(TYPE_PERM,px->max_size*sizeof(u_int),
  1500. -              new_size*sizeof(u_int));
  1501. -      }
  1502. -      px->pe = RENEW(px->pe,new_size,u_int);
  1503. -      if ( ! px->pe )
  1504. -    error(E_MEM,"px_resize");
  1505. -      px->max_size = new_size;
  1506. -   }
  1507. -   if ( px->size <= new_size )
  1508. -     /* extend permutation */
  1509. -     for ( i = px->size; i < new_size; i++ )
  1510. -       px->pe[i] = i;
  1511. -   else
  1512. -     for ( i = 0; i < new_size; i++ )
  1513. -       px->pe[i] = i;
  1514. -   
  1515. -   px->size = new_size;
  1516. -   
  1517. -   return px;
  1518. -}
  1519. -
  1520. -/* v_resize -- returns the vector x with dim new_dim
  1521. -   -- x is set to the zero vector */
  1522. -VEC    *v_resize(x,new_dim)
  1523. -VEC    *x;
  1524. -int    new_dim;
  1525. -{
  1526. -   
  1527. -   if (new_dim < 0)
  1528. -     error(E_NEG,"v_resize");
  1529. -
  1530. -   if ( ! x )
  1531. -     return v_get(new_dim);
  1532. -
  1533. -   /* nothing is changed */
  1534. -   if (new_dim == x->dim)
  1535. -     return x;
  1536. -
  1537. -   if ( x->max_dim == 0 )    /* assume that it's from sub_vec */
  1538. -     return v_get(new_dim);
  1539. -   
  1540. -   if ( new_dim > x->max_dim )
  1541. -   {
  1542. -      if (mem_info_is_on()) { 
  1543. -     mem_bytes(TYPE_VEC,x->max_dim*sizeof(Real),
  1544. -             new_dim*sizeof(Real));
  1545. -      }
  1546. -
  1547. -      x->ve = RENEW(x->ve,new_dim,Real);
  1548. -      if ( ! x->ve )
  1549. -    error(E_MEM,"v_resize");
  1550. -      x->max_dim = new_dim;
  1551. -   }
  1552. -   
  1553. -   if ( new_dim > x->dim )
  1554. -     __zero__(&(x->ve[x->dim]),new_dim - x->dim);
  1555. -   x->dim = new_dim;
  1556. -   
  1557. -   return x;
  1558. -}
  1559. -
  1560. -
  1561. -
  1562. -
  1563. -/* Varying number of arguments */
  1564. -/* other functions of this type are in sparse.c and zmemory.c */
  1565. -
  1566. -
  1567. -
  1568. -#ifdef ANSI_C
  1569. -
  1570. -
  1571. -/* To allocate memory to many arguments. 
  1572. -   The function should be called:
  1573. -   v_get_vars(dim,&x,&y,&z,...,NULL);
  1574. -   where 
  1575. -     int dim;
  1576. -     VEC *x, *y, *z,...;
  1577. -     The last argument should be NULL ! 
  1578. -     dim is the length of vectors x,y,z,...
  1579. -     returned value is equal to the number of allocated variables
  1580. -     Other gec_... functions are similar.
  1581. -*/
  1582. -
  1583. -int v_get_vars(int dim,...) 
  1584. -{
  1585. -   va_list ap;
  1586. -   int i=0;
  1587. -   VEC **par;
  1588. -   
  1589. -   va_start(ap, dim);
  1590. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  1591. -      *par = v_get(dim);
  1592. -      i++;
  1593. -   } 
  1594. -
  1595. -   va_end(ap);
  1596. -   return i;
  1597. -}
  1598. -
  1599. -
  1600. -int iv_get_vars(int dim,...) 
  1601. -{
  1602. -   va_list ap;
  1603. -   int i=0;
  1604. -   IVEC **par;
  1605. -   
  1606. -   va_start(ap, dim);
  1607. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  1608. -      *par = iv_get(dim);
  1609. -      i++;
  1610. -   } 
  1611. -
  1612. -   va_end(ap);
  1613. -   return i;
  1614. -}
  1615. -
  1616. -int m_get_vars(int m,int n,...) 
  1617. -{
  1618. -   va_list ap;
  1619. -   int i=0;
  1620. -   MAT **par;
  1621. -   
  1622. -   va_start(ap, n);
  1623. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  1624. -      *par = m_get(m,n);
  1625. -      i++;
  1626. -   } 
  1627. -
  1628. -   va_end(ap);
  1629. -   return i;
  1630. -}
  1631. -
  1632. -int px_get_vars(int dim,...) 
  1633. -{
  1634. -   va_list ap;
  1635. -   int i=0;
  1636. -   PERM **par;
  1637. -   
  1638. -   va_start(ap, dim);
  1639. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  1640. -      *par = px_get(dim);
  1641. -      i++;
  1642. -   } 
  1643. -
  1644. -   va_end(ap);
  1645. -   return i;
  1646. -}
  1647. -
  1648. -
  1649. -
  1650. -/* To resize memory for many arguments. 
  1651. -   The function should be called:
  1652. -   v_resize_vars(new_dim,&x,&y,&z,...,NULL);
  1653. -   where 
  1654. -     int new_dim;
  1655. -     VEC *x, *y, *z,...;
  1656. -     The last argument should be NULL ! 
  1657. -     rdim is the resized length of vectors x,y,z,...
  1658. -     returned value is equal to the number of allocated variables.
  1659. -     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  1660. -     argument. 
  1661. -     Other *_resize_list() functions are similar.
  1662. -*/
  1663. -
  1664. -int v_resize_vars(int new_dim,...)
  1665. -{
  1666. -   va_list ap;
  1667. -   int i=0;
  1668. -   VEC **par;
  1669. -   
  1670. -   va_start(ap, new_dim);
  1671. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  1672. -      *par = v_resize(*par,new_dim);
  1673. -      i++;
  1674. -   } 
  1675. -
  1676. -   va_end(ap);
  1677. -   return i;
  1678. -}
  1679. -
  1680. -
  1681. -
  1682. -int iv_resize_vars(int new_dim,...) 
  1683. -{
  1684. -   va_list ap;
  1685. -   int i=0;
  1686. -   IVEC **par;
  1687. -   
  1688. -   va_start(ap, new_dim);
  1689. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  1690. -      *par = iv_resize(*par,new_dim);
  1691. -      i++;
  1692. -   } 
  1693. -
  1694. -   va_end(ap);
  1695. -   return i;
  1696. -}
  1697. -
  1698. -int m_resize_vars(int m,int n,...) 
  1699. -{
  1700. -   va_list ap;
  1701. -   int i=0;
  1702. -   MAT **par;
  1703. -   
  1704. -   va_start(ap, n);
  1705. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  1706. -      *par = m_resize(*par,m,n);
  1707. -      i++;
  1708. -   } 
  1709. -
  1710. -   va_end(ap);
  1711. -   return i;
  1712. -}
  1713. -
  1714. -
  1715. -int px_resize_vars(int new_dim,...) 
  1716. -{
  1717. -   va_list ap;
  1718. -   int i=0;
  1719. -   PERM **par;
  1720. -   
  1721. -   va_start(ap, new_dim);
  1722. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  1723. -      *par = px_resize(*par,new_dim);
  1724. -      i++;
  1725. -   } 
  1726. -
  1727. -   va_end(ap);
  1728. -   return i;
  1729. -}
  1730. -
  1731. -/* To deallocate memory for many arguments. 
  1732. -   The function should be called:
  1733. -   v_free_vars(&x,&y,&z,...,NULL);
  1734. -   where 
  1735. -     VEC *x, *y, *z,...;
  1736. -     The last argument should be NULL ! 
  1737. -     There must be at least one not NULL argument.
  1738. -     returned value is equal to the number of allocated variables.
  1739. -     Returned value of x,y,z,.. is VNULL.
  1740. -     Other *_free_list() functions are similar.
  1741. -*/
  1742. -
  1743. -
  1744. -int v_free_vars(VEC **pv,...)
  1745. -{
  1746. -   va_list ap;
  1747. -   int i=1;
  1748. -   VEC **par;
  1749. -   
  1750. -   v_free(*pv);
  1751. -   *pv = VNULL;
  1752. -   va_start(ap, pv);
  1753. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  1754. -      v_free(*par); 
  1755. -      *par = VNULL;
  1756. -      i++;
  1757. -   } 
  1758. -
  1759. -   va_end(ap);
  1760. -   return i;
  1761. -}
  1762. -
  1763. -
  1764. -int iv_free_vars(IVEC **ipv,...)
  1765. -{
  1766. -   va_list ap;
  1767. -   int i=1;
  1768. -   IVEC **par;
  1769. -   
  1770. -   iv_free(*ipv);
  1771. -   *ipv = IVNULL;
  1772. -   va_start(ap, ipv);
  1773. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  1774. -      iv_free(*par); 
  1775. -      *par = IVNULL;
  1776. -      i++;
  1777. -   } 
  1778. -
  1779. -   va_end(ap);
  1780. -   return i;
  1781. -}
  1782. -
  1783. -
  1784. -int px_free_vars(PERM **vpx,...)
  1785. -{
  1786. -   va_list ap;
  1787. -   int i=1;
  1788. -   PERM **par;
  1789. -   
  1790. -   px_free(*vpx);
  1791. -   *vpx = PNULL;
  1792. -   va_start(ap, vpx);
  1793. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  1794. -      px_free(*par); 
  1795. -      *par = PNULL;
  1796. -      i++;
  1797. -   } 
  1798. -
  1799. -   va_end(ap);
  1800. -   return i;
  1801. -}
  1802. -
  1803. -int m_free_vars(MAT **va,...)
  1804. -{
  1805. -   va_list ap;
  1806. -   int i=1;
  1807. -   MAT **par;
  1808. -   
  1809. -   m_free(*va);
  1810. -   *va = MNULL;
  1811. -   va_start(ap, va);
  1812. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  1813. -      m_free(*par); 
  1814. -      *par = MNULL;
  1815. -      i++;
  1816. -   } 
  1817. -
  1818. -   va_end(ap);
  1819. -   return i;
  1820. -}
  1821. -
  1822. -
  1823. -#elif VARARGS
  1824. -/* old varargs is used */
  1825. -
  1826. -
  1827. -
  1828. -/* To allocate memory to many arguments. 
  1829. -   The function should be called:
  1830. -   v_get_vars(dim,&x,&y,&z,...,VNULL);
  1831. -   where 
  1832. -     int dim;
  1833. -     VEC *x, *y, *z,...;
  1834. -     The last argument should be VNULL ! 
  1835. -     dim is the length of vectors x,y,z,...
  1836. -*/
  1837. -
  1838. -int v_get_vars(va_alist) va_dcl
  1839. -{
  1840. -   va_list ap;
  1841. -   int dim,i=0;
  1842. -   VEC **par;
  1843. -   
  1844. -   va_start(ap);
  1845. -   dim = va_arg(ap,int);
  1846. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  1847. -      *par = v_get(dim);
  1848. -      i++;
  1849. -   } 
  1850. -
  1851. -   va_end(ap);
  1852. -   return i;
  1853. -}
  1854. -
  1855. -
  1856. -int iv_get_vars(va_alist) va_dcl
  1857. -{
  1858. -   va_list ap;
  1859. -   int i=0, dim;
  1860. -   IVEC **par;
  1861. -   
  1862. -   va_start(ap);
  1863. -   dim = va_arg(ap,int);
  1864. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  1865. -      *par = iv_get(dim);
  1866. -      i++;
  1867. -   } 
  1868. -
  1869. -   va_end(ap);
  1870. -   return i;
  1871. -}
  1872. -
  1873. -int m_get_vars(va_alist) va_dcl
  1874. -{
  1875. -   va_list ap;
  1876. -   int i=0, n, m;
  1877. -   MAT **par;
  1878. -   
  1879. -   va_start(ap);
  1880. -   m = va_arg(ap,int);
  1881. -   n = va_arg(ap,int);
  1882. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  1883. -      *par = m_get(m,n);
  1884. -      i++;
  1885. -   } 
  1886. -
  1887. -   va_end(ap);
  1888. -   return i;
  1889. -}
  1890. -
  1891. -
  1892. -
  1893. -int px_get_vars(va_alist) va_dcl
  1894. -{
  1895. -   va_list ap;
  1896. -   int i=0, dim;
  1897. -   PERM **par;
  1898. -   
  1899. -   va_start(ap);
  1900. -   dim = va_arg(ap,int);
  1901. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  1902. -      *par = px_get(dim);
  1903. -      i++;
  1904. -   } 
  1905. -
  1906. -   va_end(ap);
  1907. -   return i;
  1908. -}
  1909. -
  1910. -
  1911. -
  1912. -/* To resize memory for many arguments. 
  1913. -   The function should be called:
  1914. -   v_resize_vars(new_dim,&x,&y,&z,...,NULL);
  1915. -   where 
  1916. -     int new_dim;
  1917. -     VEC *x, *y, *z,...;
  1918. -     The last argument should be NULL ! 
  1919. -     rdim is the resized length of vectors x,y,z,...
  1920. -     returned value is equal to the number of allocated variables.
  1921. -     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  1922. -     argument. 
  1923. -     Other *_resize_list() functions are similar.
  1924. -*/
  1925. -
  1926. -int v_resize_vars(va_alist) va_dcl
  1927. -{
  1928. -   va_list ap;
  1929. -   int i=0, new_dim;
  1930. -   VEC **par;
  1931. -   
  1932. -   va_start(ap);
  1933. -   new_dim = va_arg(ap,int);
  1934. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  1935. -      *par = v_resize(*par,new_dim);
  1936. -      i++;
  1937. -   } 
  1938. -
  1939. -   va_end(ap);
  1940. -   return i;
  1941. -}
  1942. -
  1943. -
  1944. -
  1945. -int iv_resize_vars(va_alist) va_dcl
  1946. -{
  1947. -   va_list ap;
  1948. -   int i=0, new_dim;
  1949. -   IVEC **par;
  1950. -   
  1951. -   va_start(ap);
  1952. -   new_dim = va_arg(ap,int);
  1953. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  1954. -      *par = iv_resize(*par,new_dim);
  1955. -      i++;
  1956. -   } 
  1957. -
  1958. -   va_end(ap);
  1959. -   return i;
  1960. -}
  1961. -
  1962. -int m_resize_vars(va_alist) va_dcl
  1963. -{
  1964. -   va_list ap;
  1965. -   int i=0, m, n;
  1966. -   MAT **par;
  1967. -   
  1968. -   va_start(ap);
  1969. -   m = va_arg(ap,int);
  1970. -   n = va_arg(ap,int);
  1971. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  1972. -      *par = m_resize(*par,m,n);
  1973. -      i++;
  1974. -   } 
  1975. -
  1976. -   va_end(ap);
  1977. -   return i;
  1978. -}
  1979. -
  1980. -int px_resize_vars(va_alist) va_dcl
  1981. -{
  1982. -   va_list ap;
  1983. -   int i=0, new_dim;
  1984. -   PERM **par;
  1985. -   
  1986. -   va_start(ap);
  1987. -   new_dim = va_arg(ap,int);
  1988. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  1989. -      *par = px_resize(*par,new_dim);
  1990. -      i++;
  1991. -   } 
  1992. -
  1993. -   va_end(ap);
  1994. -   return i;
  1995. -}
  1996. -
  1997. -
  1998. -/* To deallocate memory for many arguments. 
  1999. -   The function should be called:
  2000. -   v_free_vars(&x,&y,&z,...,NULL);
  2001. -   where 
  2002. -     VEC *x, *y, *z,...;
  2003. -     The last argument should be NULL ! 
  2004. -     returned value is equal to the number of allocated variables.
  2005. -     Returned value of x,y,z,.. is VNULL.
  2006. -     Other *_free_list() functions are similar.
  2007. -*/
  2008. -
  2009. -
  2010. -int v_free_vars(va_alist) va_dcl
  2011. -{
  2012. -   va_list ap;
  2013. -   int i=0;
  2014. -   VEC **par;
  2015. -   
  2016. -   va_start(ap);
  2017. -   while (par = va_arg(ap,VEC **)) {   /* NULL ends the list*/
  2018. -      v_free(*par); 
  2019. -      *par = VNULL;
  2020. -      i++;
  2021. -   } 
  2022. -
  2023. -   va_end(ap);
  2024. -   return i;
  2025. -}
  2026. -
  2027. -
  2028. -
  2029. -int iv_free_vars(va_alist) va_dcl
  2030. -{
  2031. -   va_list ap;
  2032. -   int i=0;
  2033. -   IVEC **par;
  2034. -   
  2035. -   va_start(ap);
  2036. -   while (par = va_arg(ap,IVEC **)) {   /* NULL ends the list*/
  2037. -      iv_free(*par); 
  2038. -      *par = IVNULL;
  2039. -      i++;
  2040. -   } 
  2041. -
  2042. -   va_end(ap);
  2043. -   return i;
  2044. -}
  2045. -
  2046. -
  2047. -int px_free_vars(va_alist) va_dcl
  2048. -{
  2049. -   va_list ap;
  2050. -   int i=0;
  2051. -   PERM **par;
  2052. -   
  2053. -   va_start(ap);
  2054. -   while (par = va_arg(ap,PERM **)) {   /* NULL ends the list*/
  2055. -      px_free(*par); 
  2056. -      *par = PNULL;
  2057. -      i++;
  2058. -   } 
  2059. -
  2060. -   va_end(ap);
  2061. -   return i;
  2062. -}
  2063. -
  2064. -int m_free_vars(va_alist) va_dcl
  2065. -{
  2066. -   va_list ap;
  2067. -   int i=0;
  2068. -   MAT **par;
  2069. -   
  2070. -   va_start(ap);
  2071. -   while (par = va_arg(ap,MAT **)) {   /* NULL ends the list*/
  2072. -      m_free(*par); 
  2073. -      *par = MNULL;
  2074. -      i++;
  2075. -   } 
  2076. -
  2077. -   va_end(ap);
  2078. -   return i;
  2079. -}
  2080. -
  2081. -
  2082. -
  2083. -#endif /* VARARGS */
  2084. -  
  2085. -
  2086. //GO.SYSIN DD memory.c
  2087. echo vecop.c 1>&2
  2088. sed >vecop.c <<'//GO.SYSIN DD vecop.c' 's/^-//'
  2089. -
  2090. -/**************************************************************************
  2091. -**
  2092. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2093. -**
  2094. -**                 Meschach Library
  2095. -** 
  2096. -** This Meschach Library is provided "as is" without any express 
  2097. -** or implied warranty of any kind with respect to this software. 
  2098. -** In particular the authors shall not be liable for any direct, 
  2099. -** indirect, special, incidental or consequential damages arising 
  2100. -** in any way from use of the software.
  2101. -** 
  2102. -** Everyone is granted permission to copy, modify and redistribute this
  2103. -** Meschach Library, provided:
  2104. -**  1.  All copies contain this copyright notice.
  2105. -**  2.  All modified copies shall carry a notice stating who
  2106. -**      made the last modification and the date of such modification.
  2107. -**  3.  No charge is made for this software or works derived from it.  
  2108. -**      This clause shall not be construed as constraining other software
  2109. -**      distributed on the same medium as this software, nor is a
  2110. -**      distribution fee considered a charge.
  2111. -**
  2112. -***************************************************************************/
  2113. -
  2114. -
  2115. -/* vecop.c 1.3 8/18/87 */
  2116. -
  2117. -#include    <stdio.h>
  2118. -#include    "matrix.h"
  2119. -
  2120. -static    char    rcsid[] = "$Id: vecop.c,v 1.4 1994/03/08 05:50:39 des Exp $";
  2121. -
  2122. -
  2123. -/* _in_prod -- inner product of two vectors from i0 downwards */
  2124. -double    _in_prod(a,b,i0)
  2125. -VEC    *a,*b;
  2126. -u_int    i0;
  2127. -{
  2128. -    u_int    limit;
  2129. -    /* Real    *a_v, *b_v; */
  2130. -    /* register Real    sum; */
  2131. -
  2132. -    if ( a==(VEC *)NULL || b==(VEC *)NULL )
  2133. -        error(E_NULL,"_in_prod");
  2134. -    limit = min(a->dim,b->dim);
  2135. -    if ( i0 > limit )
  2136. -        error(E_BOUNDS,"_in_prod");
  2137. -
  2138. -    return __ip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0));
  2139. -    /*****************************************
  2140. -    a_v = &(a->ve[i0]);        b_v = &(b->ve[i0]);
  2141. -    for ( i=i0; i<limit; i++ )
  2142. -        sum += a_v[i]*b_v[i];
  2143. -        sum += (*a_v++)*(*b_v++);
  2144. -
  2145. -    return (double)sum;
  2146. -    ******************************************/
  2147. -}
  2148. -
  2149. -/* sv_mlt -- scalar-vector multiply -- may be in-situ */
  2150. -VEC    *sv_mlt(scalar,vector,out)
  2151. -double    scalar;
  2152. -VEC    *vector,*out;
  2153. -{
  2154. -    /* u_int    dim, i; */
  2155. -    /* Real    *out_ve, *vec_ve; */
  2156. -
  2157. -    if ( vector==(VEC *)NULL )
  2158. -        error(E_NULL,"sv_mlt");
  2159. -    if ( out==(VEC *)NULL || out->dim != vector->dim )
  2160. -        out = v_resize(out,vector->dim);
  2161. -    if ( scalar == 0.0 )
  2162. -        return v_zero(out);
  2163. -    if ( scalar == 1.0 )
  2164. -        return v_copy(vector,out);
  2165. -
  2166. -    __smlt__(vector->ve,(double)scalar,out->ve,(int)(vector->dim));
  2167. -    /**************************************************
  2168. -    dim = vector->dim;
  2169. -    out_ve = out->ve;    vec_ve = vector->ve;
  2170. -    for ( i=0; i<dim; i++ )
  2171. -        out->ve[i] = scalar*vector->ve[i];
  2172. -        (*out_ve++) = scalar*(*vec_ve++);
  2173. -    **************************************************/
  2174. -    return (out);
  2175. -}
  2176. -
  2177. -/* v_add -- vector addition -- may be in-situ */
  2178. -VEC    *v_add(vec1,vec2,out)
  2179. -VEC    *vec1,*vec2,*out;
  2180. -{
  2181. -    u_int    dim;
  2182. -    /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  2183. -
  2184. -    if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  2185. -        error(E_NULL,"v_add");
  2186. -    if ( vec1->dim != vec2->dim )
  2187. -        error(E_SIZES,"v_add");
  2188. -    if ( out==(VEC *)NULL || out->dim != vec1->dim )
  2189. -        out = v_resize(out,vec1->dim);
  2190. -    dim = vec1->dim;
  2191. -    __add__(vec1->ve,vec2->ve,out->ve,(int)dim);
  2192. -    /************************************************************
  2193. -    out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  2194. -    for ( i=0; i<dim; i++ )
  2195. -        out->ve[i] = vec1->ve[i]+vec2->ve[i];
  2196. -        (*out_ve++) = (*vec1_ve++) + (*vec2_ve++);
  2197. -    ************************************************************/
  2198. -
  2199. -    return (out);
  2200. -}
  2201. -
  2202. -/* v_mltadd -- scalar/vector multiplication and addition
  2203. -        -- out = v1 + scale.v2        */
  2204. -VEC    *v_mltadd(v1,v2,scale,out)
  2205. -VEC    *v1,*v2,*out;
  2206. -double    scale;
  2207. -{
  2208. -    /* register u_int    dim, i; */
  2209. -    /* Real    *out_ve, *v1_ve, *v2_ve; */
  2210. -
  2211. -    if ( v1==(VEC *)NULL || v2==(VEC *)NULL )
  2212. -        error(E_NULL,"v_mltadd");
  2213. -    if ( v1->dim != v2->dim )
  2214. -        error(E_SIZES,"v_mltadd");
  2215. -    if ( scale == 0.0 )
  2216. -        return v_copy(v1,out);
  2217. -    if ( scale == 1.0 )
  2218. -        return v_add(v1,v2,out);
  2219. -
  2220. -    if ( v2 != out )
  2221. -    {
  2222. -        tracecatch(out = v_copy(v1,out),"v_mltadd");
  2223. -
  2224. -        /* dim = v1->dim; */
  2225. -        __mltadd__(out->ve,v2->ve,scale,(int)(v1->dim));
  2226. -    }
  2227. -    else
  2228. -    {
  2229. -        tracecatch(out = sv_mlt(scale,v2,out),"v_mltadd");
  2230. -        out = v_add(v1,out,out);
  2231. -    }
  2232. -    /************************************************************
  2233. -    out_ve = out->ve;    v1_ve = v1->ve;        v2_ve = v2->ve;
  2234. -    for ( i=0; i < dim ; i++ )
  2235. -        out->ve[i] = v1->ve[i] + scale*v2->ve[i];
  2236. -        (*out_ve++) = (*v1_ve++) + scale*(*v2_ve++);
  2237. -    ************************************************************/
  2238. -
  2239. -    return (out);
  2240. -}
  2241. -
  2242. -/* v_sub -- vector subtraction -- may be in-situ */
  2243. -VEC    *v_sub(vec1,vec2,out)
  2244. -VEC    *vec1,*vec2,*out;
  2245. -{
  2246. -    /* u_int    i, dim; */
  2247. -    /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  2248. -
  2249. -    if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  2250. -        error(E_NULL,"v_sub");
  2251. -    if ( vec1->dim != vec2->dim )
  2252. -        error(E_SIZES,"v_sub");
  2253. -    if ( out==(VEC *)NULL || out->dim != vec1->dim )
  2254. -        out = v_resize(out,vec1->dim);
  2255. -
  2256. -    __sub__(vec1->ve,vec2->ve,out->ve,(int)(vec1->dim));
  2257. -    /************************************************************
  2258. -    dim = vec1->dim;
  2259. -    out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  2260. -    for ( i=0; i<dim; i++ )
  2261. -        out->ve[i] = vec1->ve[i]-vec2->ve[i];
  2262. -        (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
  2263. -    ************************************************************/
  2264. -
  2265. -    return (out);
  2266. -}
  2267. -
  2268. -/* v_map -- maps function f over components of x: out[i] = f(x[i])
  2269. -    -- _v_map sets out[i] = f(params,x[i]) */
  2270. -VEC    *v_map(f,x,out)
  2271. -#ifdef PROTOTYPES_IN_STRUCT
  2272. -double    (*f)(double);
  2273. -#else
  2274. -double    (*f)();
  2275. -#endif
  2276. -VEC    *x, *out;
  2277. -{
  2278. -    Real    *x_ve, *out_ve;
  2279. -    int    i, dim;
  2280. -
  2281. -    if ( ! x || ! f )
  2282. -        error(E_NULL,"v_map");
  2283. -    if ( ! out || out->dim != x->dim )
  2284. -        out = v_resize(out,x->dim);
  2285. -
  2286. -    dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  2287. -    for ( i = 0; i < dim; i++ )
  2288. -        *out_ve++ = (*f)(*x_ve++);
  2289. -
  2290. -    return out;
  2291. -}
  2292. -
  2293. -VEC    *_v_map(f,params,x,out)
  2294. -#ifdef PROTOTYPES_IN_STRUCT
  2295. -double    (*f)(void *,double);
  2296. -#else
  2297. -double    (*f)();
  2298. -#endif
  2299. -VEC    *x, *out;
  2300. -void    *params;
  2301. -{
  2302. -    Real    *x_ve, *out_ve;
  2303. -    int    i, dim;
  2304. -
  2305. -    if ( ! x || ! f )
  2306. -        error(E_NULL,"_v_map");
  2307. -    if ( ! out || out->dim != x->dim )
  2308. -        out = v_resize(out,x->dim);
  2309. -
  2310. -    dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  2311. -    for ( i = 0; i < dim; i++ )
  2312. -        *out_ve++ = (*f)(params,*x_ve++);
  2313. -
  2314. -    return out;
  2315. -}
  2316. -
  2317. -/* v_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */
  2318. -VEC    *v_lincomb(n,v,a,out)
  2319. -int    n;    /* number of a's and v's */
  2320. -Real    a[];
  2321. -VEC    *v[], *out;
  2322. -{
  2323. -    int    i;
  2324. -
  2325. -    if ( ! a || ! v )
  2326. -        error(E_NULL,"v_lincomb");
  2327. -    if ( n <= 0 )
  2328. -        return VNULL;
  2329. -
  2330. -    for ( i = 1; i < n; i++ )
  2331. -        if ( out == v[i] )
  2332. -            error(E_INSITU,"v_lincomb");
  2333. -
  2334. -    out = sv_mlt(a[0],v[0],out);
  2335. -    for ( i = 1; i < n; i++ )
  2336. -    {
  2337. -        if ( ! v[i] )
  2338. -            error(E_NULL,"v_lincomb");
  2339. -        if ( v[i]->dim != out->dim )
  2340. -            error(E_SIZES,"v_lincomb");
  2341. -        out = v_mltadd(out,v[i],a[i],out);
  2342. -    }
  2343. -
  2344. -    return out;
  2345. -}
  2346. -
  2347. -
  2348. -
  2349. -#ifdef ANSI_C
  2350. -
  2351. -/* v_linlist -- linear combinations taken from a list of arguments;
  2352. -   calling:
  2353. -      v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  2354. -   where vi are vectors (VEC *) and ai are numbers (double)
  2355. -*/
  2356. -VEC  *v_linlist(VEC *out,VEC *v1,double a1,...)
  2357. -{
  2358. -   va_list ap;
  2359. -   VEC *par;
  2360. -   double a_par;
  2361. -
  2362. -   if ( ! v1 )
  2363. -     return VNULL;
  2364. -   
  2365. -   va_start(ap, a1);
  2366. -   out = sv_mlt(a1,v1,out);
  2367. -   
  2368. -   while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  2369. -      a_par = va_arg(ap,double);
  2370. -      if (a_par == 0.0) continue;
  2371. -      if ( out == par )        
  2372. -    error(E_INSITU,"v_linlist");
  2373. -      if ( out->dim != par->dim )    
  2374. -    error(E_SIZES,"v_linlist");
  2375. -
  2376. -      if (a_par == 1.0)
  2377. -    out = v_add(out,par,out);
  2378. -      else if (a_par == -1.0)
  2379. -    out = v_sub(out,par,out);
  2380. -      else
  2381. -    out = v_mltadd(out,par,a_par,out); 
  2382. -   } 
  2383. -   
  2384. -   va_end(ap);
  2385. -   return out;
  2386. -}
  2387. - 
  2388. -#elif VARARGS
  2389. -
  2390. -
  2391. -/* v_linlist -- linear combinations taken from a list of arguments;
  2392. -   calling:
  2393. -      v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  2394. -   where vi are vectors (VEC *) and ai are numbers (double)
  2395. -*/
  2396. -VEC  *v_linlist(va_alist) va_dcl
  2397. -{
  2398. -   va_list ap;
  2399. -   VEC *par, *out;
  2400. -   double a_par;
  2401. -
  2402. -   va_start(ap);
  2403. -   out = va_arg(ap,VEC *);
  2404. -   par = va_arg(ap,VEC *);
  2405. -   if ( ! par ) {
  2406. -      va_end(ap);
  2407. -      return VNULL;
  2408. -   }
  2409. -   
  2410. -   a_par = va_arg(ap,double);
  2411. -   out = sv_mlt(a_par,par,out);
  2412. -   
  2413. -   while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  2414. -      a_par = va_arg(ap,double);
  2415. -      if (a_par == 0.0) continue;
  2416. -      if ( out == par )        
  2417. -    error(E_INSITU,"v_linlist");
  2418. -      if ( out->dim != par->dim )    
  2419. -    error(E_SIZES,"v_linlist");
  2420. -
  2421. -      if (a_par == 1.0)
  2422. -    out = v_add(out,par,out);
  2423. -      else if (a_par == -1.0)
  2424. -    out = v_sub(out,par,out);
  2425. -      else
  2426. -    out = v_mltadd(out,par,a_par,out); 
  2427. -   } 
  2428. -   
  2429. -   va_end(ap);
  2430. -   return out;
  2431. -}
  2432. -
  2433. -#endif
  2434. -  
  2435. -
  2436. -
  2437. -
  2438. -
  2439. -/* v_star -- computes componentwise (Hadamard) product of x1 and x2
  2440. -    -- result out is returned */
  2441. -VEC    *v_star(x1, x2, out)
  2442. -VEC    *x1, *x2, *out;
  2443. -{
  2444. -    int        i;
  2445. -
  2446. -    if ( ! x1 || ! x2 )
  2447. -    error(E_NULL,"v_star");
  2448. -    if ( x1->dim != x2->dim )
  2449. -    error(E_SIZES,"v_star");
  2450. -    out = v_resize(out,x1->dim);
  2451. -
  2452. -    for ( i = 0; i < x1->dim; i++ )
  2453. -    out->ve[i] = x1->ve[i] * x2->ve[i];
  2454. -
  2455. -    return out;
  2456. -}
  2457. -
  2458. -/* v_slash -- computes componentwise ratio of x2 and x1
  2459. -    -- out[i] = x2[i] / x1[i]
  2460. -    -- if x1[i] == 0 for some i, then raise E_SING error
  2461. -    -- result out is returned */
  2462. -VEC    *v_slash(x1, x2, out)
  2463. -VEC    *x1, *x2, *out;
  2464. -{
  2465. -    int        i;
  2466. -    Real    tmp;
  2467. -
  2468. -    if ( ! x1 || ! x2 )
  2469. -    error(E_NULL,"v_slash");
  2470. -    if ( x1->dim != x2->dim )
  2471. -    error(E_SIZES,"v_slash");
  2472. -    out = v_resize(out,x1->dim);
  2473. -
  2474. -    for ( i = 0; i < x1->dim; i++ )
  2475. -    {
  2476. -    tmp = x1->ve[i];
  2477. -    if ( tmp == 0.0 )
  2478. -        error(E_SING,"v_slash");
  2479. -    out->ve[i] = x2->ve[i] / tmp;
  2480. -    }
  2481. -
  2482. -    return out;
  2483. -}
  2484. -
  2485. -/* v_min -- computes minimum component of x, which is returned
  2486. -    -- also sets min_idx to the index of this minimum */
  2487. -double    v_min(x, min_idx)
  2488. -VEC    *x;
  2489. -int    *min_idx;
  2490. -{
  2491. -    int        i, i_min;
  2492. -    Real    min_val, tmp;
  2493. -
  2494. -    if ( ! x )
  2495. -    error(E_NULL,"v_min");
  2496. -    if ( x->dim <= 0 )
  2497. -    error(E_SIZES,"v_min");
  2498. -    i_min = 0;
  2499. -    min_val = x->ve[0];
  2500. -    for ( i = 1; i < x->dim; i++ )
  2501. -    {
  2502. -    tmp = x->ve[i];
  2503. -    if ( tmp < min_val )
  2504. -    {
  2505. -        min_val = tmp;
  2506. -        i_min = i;
  2507. -    }
  2508. -    }
  2509. -
  2510. -    if ( min_idx != NULL )
  2511. -    *min_idx = i_min;
  2512. -    return min_val;
  2513. -}
  2514. -
  2515. -/* v_max -- computes maximum component of x, which is returned
  2516. -    -- also sets max_idx to the index of this maximum */
  2517. -double    v_max(x, max_idx)
  2518. -VEC    *x;
  2519. -int    *max_idx;
  2520. -{
  2521. -    int        i, i_max;
  2522. -    Real    max_val, tmp;
  2523. -
  2524. -    if ( ! x )
  2525. -    error(E_NULL,"v_max");
  2526. -    if ( x->dim <= 0 )
  2527. -    error(E_SIZES,"v_max");
  2528. -    i_max = 0;
  2529. -    max_val = x->ve[0];
  2530. -    for ( i = 1; i < x->dim; i++ )
  2531. -    {
  2532. -    tmp = x->ve[i];
  2533. -    if ( tmp > max_val )
  2534. -    {
  2535. -        max_val = tmp;
  2536. -        i_max = i;
  2537. -    }
  2538. -    }
  2539. -
  2540. -    if ( max_idx != NULL )
  2541. -    *max_idx = i_max;
  2542. -    return max_val;
  2543. -}
  2544. -
  2545. -#define    MAX_STACK    60
  2546. -
  2547. -
  2548. -/* v_sort -- sorts vector x, and generates permutation that gives the order
  2549. -    of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
  2550. -    the permutation is order = [2, 0, 1].
  2551. -    -- if order is NULL on entry then it is ignored
  2552. -    -- the sorted vector x is returned */
  2553. -VEC    *v_sort(x, order)
  2554. -VEC    *x;
  2555. -PERM    *order;
  2556. -{
  2557. -    Real    *x_ve, tmp, v;
  2558. -    /* int        *order_pe; */
  2559. -    int        dim, i, j, l, r, tmp_i;
  2560. -    int        stack[MAX_STACK], sp;
  2561. -
  2562. -    if ( ! x )
  2563. -    error(E_NULL,"v_sort");
  2564. -    if ( order != PNULL && order->size != x->dim )
  2565. -    order = px_resize(order, x->dim);
  2566. -
  2567. -    x_ve = x->ve;
  2568. -    dim = x->dim;
  2569. -    if ( order != PNULL )
  2570. -    px_ident(order);
  2571. -
  2572. -    if ( dim <= 1 )
  2573. -    return x;
  2574. -
  2575. -    /* using quicksort algorithm in Sedgewick,
  2576. -       "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
  2577. -    sp = 0;
  2578. -    l = 0;    r = dim-1;    v = x_ve[0];
  2579. -    for ( ; ; )
  2580. -    {
  2581. -    while ( r > l )
  2582. -    {
  2583. -        /* "i = partition(x_ve,l,r);" */
  2584. -        v = x_ve[r];
  2585. -        i = l-1;
  2586. -        j = r;
  2587. -        for ( ; ; )
  2588. -        {
  2589. -        while ( x_ve[++i] < v )
  2590. -            ;
  2591. -        while ( x_ve[--j] > v )
  2592. -            ;
  2593. -        if ( i >= j )    break;
  2594. -        
  2595. -        tmp = x_ve[i];
  2596. -        x_ve[i] = x_ve[j];
  2597. -        x_ve[j] = tmp;
  2598. -        if ( order != PNULL )
  2599. -        {
  2600. -            tmp_i = order->pe[i];
  2601. -            order->pe[i] = order->pe[j];
  2602. -            order->pe[j] = tmp_i;
  2603. -        }
  2604. -        }
  2605. -        tmp = x_ve[i];
  2606. -        x_ve[i] = x_ve[r];
  2607. -        x_ve[r] = tmp;
  2608. -        if ( order != PNULL )
  2609. -        {
  2610. -        tmp_i = order->pe[i];
  2611. -        order->pe[i] = order->pe[r];
  2612. -        order->pe[r] = tmp_i;
  2613. -        }
  2614. -
  2615. -        if ( i-l > r-i )
  2616. -        {   stack[sp++] = l;   stack[sp++] = i-1;   l = i+1;   }
  2617. -        else
  2618. -        {   stack[sp++] = i+1;   stack[sp++] = r;   r = i-1;   }
  2619. -    }
  2620. -
  2621. -    /* recursion elimination */
  2622. -    if ( sp == 0 )
  2623. -        break;
  2624. -    r = stack[--sp];
  2625. -    l = stack[--sp];
  2626. -    }
  2627. -
  2628. -    return x;
  2629. -}
  2630. -
  2631. -/* v_sum -- returns sum of entries of a vector */
  2632. -double    v_sum(x)
  2633. -VEC    *x;
  2634. -{
  2635. -    int        i;
  2636. -    Real    sum;
  2637. -
  2638. -    if ( ! x )
  2639. -    error(E_NULL,"v_sum");
  2640. -
  2641. -    sum = 0.0;
  2642. -    for ( i = 0; i < x->dim; i++ )
  2643. -    sum += x->ve[i];
  2644. -
  2645. -    return sum;
  2646. -}
  2647. -
  2648. -/* v_conv -- computes convolution product of two vectors */
  2649. -VEC    *v_conv(x1, x2, out)
  2650. -VEC    *x1, *x2, *out;
  2651. -{
  2652. -    int        i;
  2653. -
  2654. -    if ( ! x1 || ! x2 )
  2655. -    error(E_NULL,"v_conv");
  2656. -    if ( x1 == out || x2 == out )
  2657. -    error(E_INSITU,"v_conv");
  2658. -    if ( x1->dim == 0 || x2->dim == 0 )
  2659. -    return out = v_resize(out,0);
  2660. -
  2661. -    out = v_resize(out,x1->dim + x2->dim - 1);
  2662. -    v_zero(out);
  2663. -    for ( i = 0; i < x1->dim; i++ )
  2664. -    __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim);
  2665. -
  2666. -    return out;
  2667. -}
  2668. -
  2669. -/* v_pconv -- computes a periodic convolution product
  2670. -    -- the period is the dimension of x2 */
  2671. -VEC    *v_pconv(x1, x2, out)
  2672. -VEC    *x1, *x2, *out;
  2673. -{
  2674. -    int        i;
  2675. -
  2676. -    if ( ! x1 || ! x2 )
  2677. -    error(E_NULL,"v_pconv");
  2678. -    if ( x1 == out || x2 == out )
  2679. -    error(E_INSITU,"v_pconv");
  2680. -    out = v_resize(out,x2->dim);
  2681. -    if ( x2->dim == 0 )
  2682. -    return out;
  2683. -
  2684. -    v_zero(out);
  2685. -    for ( i = 0; i < x1->dim; i++ )
  2686. -    {
  2687. -    __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim - i);
  2688. -    if ( i > 0 )
  2689. -        __mltadd__(out->ve,&(x2->ve[x2->dim - i]),x1->ve[i],i);
  2690. -    }
  2691. -
  2692. -    return out;
  2693. -}
  2694. //GO.SYSIN DD vecop.c
  2695. echo matop.c 1>&2
  2696. sed >matop.c <<'//GO.SYSIN DD matop.c' 's/^-//'
  2697. -
  2698. -/**************************************************************************
  2699. -**
  2700. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2701. -**
  2702. -**                 Meschach Library
  2703. -** 
  2704. -** This Meschach Library is provided "as is" without any express 
  2705. -** or implied warranty of any kind with respect to this software. 
  2706. -** In particular the authors shall not be liable for any direct, 
  2707. -** indirect, special, incidental or consequential damages arising 
  2708. -** in any way from use of the software.
  2709. -** 
  2710. -** Everyone is granted permission to copy, modify and redistribute this
  2711. -** Meschach Library, provided:
  2712. -**  1.  All copies contain this copyright notice.
  2713. -**  2.  All modified copies shall carry a notice stating who
  2714. -**      made the last modification and the date of such modification.
  2715. -**  3.  No charge is made for this software or works derived from it.  
  2716. -**      This clause shall not be construed as constraining other software
  2717. -**      distributed on the same medium as this software, nor is a
  2718. -**      distribution fee considered a charge.
  2719. -**
  2720. -***************************************************************************/
  2721. -
  2722. -
  2723. -/* matop.c 1.3 11/25/87 */
  2724. -
  2725. -
  2726. -#include    <stdio.h>
  2727. -#include    "matrix.h"
  2728. -
  2729. -static    char    rcsid[] = "$Id: matop.c,v 1.3 1994/01/13 05:30:28 des Exp $";
  2730. -
  2731. -
  2732. -/* m_add -- matrix addition -- may be in-situ */
  2733. -MAT    *m_add(mat1,mat2,out)
  2734. -MAT    *mat1,*mat2,*out;
  2735. -{
  2736. -    u_int    m,n,i;
  2737. -
  2738. -    if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  2739. -        error(E_NULL,"m_add");
  2740. -    if ( mat1->m != mat2->m || mat1->n != mat2->n )
  2741. -        error(E_SIZES,"m_add");
  2742. -    if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  2743. -        out = m_resize(out,mat1->m,mat1->n);
  2744. -    m = mat1->m;    n = mat1->n;
  2745. -    for ( i=0; i<m; i++ )
  2746. -    {
  2747. -        __add__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  2748. -        /**************************************************
  2749. -        for ( j=0; j<n; j++ )
  2750. -            out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  2751. -        **************************************************/
  2752. -    }
  2753. -
  2754. -    return (out);
  2755. -}
  2756. -
  2757. -/* m_sub -- matrix subtraction -- may be in-situ */
  2758. -MAT    *m_sub(mat1,mat2,out)
  2759. -MAT    *mat1,*mat2,*out;
  2760. -{
  2761. -    u_int    m,n,i;
  2762. -
  2763. -    if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  2764. -        error(E_NULL,"m_sub");
  2765. -    if ( mat1->m != mat2->m || mat1->n != mat2->n )
  2766. -        error(E_SIZES,"m_sub");
  2767. -    if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  2768. -        out = m_resize(out,mat1->m,mat1->n);
  2769. -    m = mat1->m;    n = mat1->n;
  2770. -    for ( i=0; i<m; i++ )
  2771. -    {
  2772. -        __sub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  2773. -        /**************************************************
  2774. -        for ( j=0; j<n; j++ )
  2775. -            out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  2776. -        **************************************************/
  2777. -    }
  2778. -
  2779. -    return (out);
  2780. -}
  2781. -
  2782. -/* m_mlt -- matrix-matrix multiplication */
  2783. -MAT    *m_mlt(A,B,OUT)
  2784. -MAT    *A,*B,*OUT;
  2785. -{
  2786. -    u_int    i, /* j, */ k, m, n, p;
  2787. -    Real    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  2788. -
  2789. -    if ( A==(MAT *)NULL || B==(MAT *)NULL )
  2790. -        error(E_NULL,"m_mlt");
  2791. -    if ( A->n != B->m )
  2792. -        error(E_SIZES,"m_mlt");
  2793. -    if ( A == OUT || B == OUT )
  2794. -        error(E_INSITU,"m_mlt");
  2795. -    m = A->m;    n = A->n;    p = B->n;
  2796. -    A_v = A->me;        B_v = B->me;
  2797. -
  2798. -    if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n )
  2799. -        OUT = m_resize(OUT,A->m,B->n);
  2800. -
  2801. -/****************************************************************
  2802. -    for ( i=0; i<m; i++ )
  2803. -        for  ( j=0; j<p; j++ )
  2804. -        {
  2805. -            sum = 0.0;
  2806. -            for ( k=0; k<n; k++ )
  2807. -                sum += A_v[i][k]*B_v[k][j];
  2808. -            OUT->me[i][j] = sum;
  2809. -        }
  2810. -****************************************************************/
  2811. -    m_zero(OUT);
  2812. -    for ( i=0; i<m; i++ )
  2813. -        for ( k=0; k<n; k++ )
  2814. -        {
  2815. -            if ( A_v[i][k] != 0.0 )
  2816. -                __mltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p);
  2817. -            /**************************************************
  2818. -            B_row = B_v[k];    OUT_row = OUT->me[i];
  2819. -            for ( j=0; j<p; j++ )
  2820. -            (*OUT_row++) += tmp*(*B_row++);
  2821. -            **************************************************/
  2822. -        }
  2823. -
  2824. -    return OUT;
  2825. -}
  2826. -
  2827. -/* mmtr_mlt -- matrix-matrix transposed multiplication
  2828. -    -- A.B^T is returned, and stored in OUT */
  2829. -MAT    *mmtr_mlt(A,B,OUT)
  2830. -MAT    *A, *B, *OUT;
  2831. -{
  2832. -    int    i, j, limit;
  2833. -    /* Real    *A_row, *B_row, sum; */
  2834. -
  2835. -    if ( ! A || ! B )
  2836. -        error(E_NULL,"mmtr_mlt");
  2837. -    if ( A == OUT || B == OUT )
  2838. -        error(E_INSITU,"mmtr_mlt");
  2839. -    if ( A->n != B->n )
  2840. -        error(E_SIZES,"mmtr_mlt");
  2841. -    if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  2842. -        OUT = m_resize(OUT,A->m,B->m);
  2843. -
  2844. -    limit = A->n;
  2845. -    for ( i = 0; i < A->m; i++ )
  2846. -        for ( j = 0; j < B->m; j++ )
  2847. -        {
  2848. -            OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit);
  2849. -            /**************************************************
  2850. -            sum = 0.0;
  2851. -            A_row = A->me[i];
  2852. -            B_row = B->me[j];
  2853. -            for ( k = 0; k < limit; k++ )
  2854. -            sum += (*A_row++)*(*B_row++);
  2855. -            OUT->me[i][j] = sum;
  2856. -            **************************************************/
  2857. -        }
  2858. -
  2859. -    return OUT;
  2860. -}
  2861. -
  2862. -/* mtrm_mlt -- matrix transposed-matrix multiplication
  2863. -    -- A^T.B is returned, result stored in OUT */
  2864. -MAT    *mtrm_mlt(A,B,OUT)
  2865. -MAT    *A, *B, *OUT;
  2866. -{
  2867. -    int    i, k, limit;
  2868. -    /* Real    *B_row, *OUT_row, multiplier; */
  2869. -
  2870. -    if ( ! A || ! B )
  2871. -        error(E_NULL,"mmtr_mlt");
  2872. -    if ( A == OUT || B == OUT )
  2873. -        error(E_INSITU,"mtrm_mlt");
  2874. -    if ( A->m != B->m )
  2875. -        error(E_SIZES,"mmtr_mlt");
  2876. -    if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  2877. -        OUT = m_resize(OUT,A->n,B->n);
  2878. -
  2879. -    limit = B->n;
  2880. -    m_zero(OUT);
  2881. -    for ( k = 0; k < A->m; k++ )
  2882. -        for ( i = 0; i < A->n; i++ )
  2883. -        {
  2884. -            if ( A->me[k][i] != 0.0 )
  2885. -            __mltadd__(OUT->me[i],B->me[k],A->me[k][i],(int)limit);
  2886. -            /**************************************************
  2887. -            multiplier = A->me[k][i];
  2888. -            OUT_row = OUT->me[i];
  2889. -            B_row   = B->me[k];
  2890. -            for ( j = 0; j < limit; j++ )
  2891. -            *(OUT_row++) += multiplier*(*B_row++);
  2892. -            **************************************************/
  2893. -        }
  2894. -
  2895. -    return OUT;
  2896. -}
  2897. -
  2898. -/* mv_mlt -- matrix-vector multiplication 
  2899. -        -- Note: b is treated as a column vector */
  2900. -VEC    *mv_mlt(A,b,out)
  2901. -MAT    *A;
  2902. -VEC    *b,*out;
  2903. -{
  2904. -    u_int    i, m, n;
  2905. -    Real    **A_v, *b_v /*, *A_row */;
  2906. -    /* register Real    sum; */
  2907. -
  2908. -    if ( A==(MAT *)NULL || b==(VEC *)NULL )
  2909. -        error(E_NULL,"mv_mlt");
  2910. -    if ( A->n != b->dim )
  2911. -        error(E_SIZES,"mv_mlt");
  2912. -    if ( b == out )
  2913. -        error(E_INSITU,"mv_mlt");
  2914. -    if ( out == (VEC *)NULL || out->dim != A->m )
  2915. -        out = v_resize(out,A->m);
  2916. -
  2917. -    m = A->m;        n = A->n;
  2918. -    A_v = A->me;        b_v = b->ve;
  2919. -    for ( i=0; i<m; i++ )
  2920. -    {
  2921. -        /* for ( j=0; j<n; j++ )
  2922. -            sum += A_v[i][j]*b_v[j]; */
  2923. -        out->ve[i] = __ip__(A_v[i],b_v,(int)n);
  2924. -        /**************************************************
  2925. -        A_row = A_v[i];        b_v = b->ve;
  2926. -        for ( j=0; j<n; j++ )
  2927. -            sum += (*A_row++)*(*b_v++);
  2928. -        out->ve[i] = sum;
  2929. -        **************************************************/
  2930. -    }
  2931. -
  2932. -    return out;
  2933. -}
  2934. -
  2935. -/* sm_mlt -- scalar-matrix multiply -- may be in-situ */
  2936. -MAT    *sm_mlt(scalar,matrix,out)
  2937. -double    scalar;
  2938. -MAT    *matrix,*out;
  2939. -{
  2940. -    u_int    m,n,i;
  2941. -
  2942. -    if ( matrix==(MAT *)NULL )
  2943. -        error(E_NULL,"sm_mlt");
  2944. -    if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n )
  2945. -        out = m_resize(out,matrix->m,matrix->n);
  2946. -    m = matrix->m;    n = matrix->n;
  2947. -    for ( i=0; i<m; i++ )
  2948. -        __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
  2949. -        /**************************************************
  2950. -        for ( j=0; j<n; j++ )
  2951. -            out->me[i][j] = scalar*matrix->me[i][j];
  2952. -        **************************************************/
  2953. -    return (out);
  2954. -}
  2955. -
  2956. -/* vm_mlt -- vector-matrix multiplication 
  2957. -        -- Note: b is treated as a row vector */
  2958. -VEC    *vm_mlt(A,b,out)
  2959. -MAT    *A;
  2960. -VEC    *b,*out;
  2961. -{
  2962. -    u_int    j,m,n;
  2963. -    /* Real    sum,**A_v,*b_v; */
  2964. -
  2965. -    if ( A==(MAT *)NULL || b==(VEC *)NULL )
  2966. -        error(E_NULL,"vm_mlt");
  2967. -    if ( A->m != b->dim )
  2968. -        error(E_SIZES,"vm_mlt");
  2969. -    if ( b == out )
  2970. -        error(E_INSITU,"vm_mlt");
  2971. -    if ( out == (VEC *)NULL || out->dim != A->n )
  2972. -        out = v_resize(out,A->n);
  2973. -
  2974. -    m = A->m;        n = A->n;
  2975. -
  2976. -    v_zero(out);
  2977. -    for ( j = 0; j < m; j++ )
  2978. -        if ( b->ve[j] != 0.0 )
  2979. -            __mltadd__(out->ve,A->me[j],b->ve[j],(int)n);
  2980. -    /**************************************************
  2981. -    A_v = A->me;        b_v = b->ve;
  2982. -    for ( j=0; j<n; j++ )
  2983. -    {
  2984. -        sum = 0.0;
  2985. -        for ( i=0; i<m; i++ )
  2986. -            sum += b_v[i]*A_v[i][j];
  2987. -        out->ve[j] = sum;
  2988. -    }
  2989. -    **************************************************/
  2990. -
  2991. -    return out;
  2992. -}
  2993. -
  2994. -/* m_transp -- transpose matrix */
  2995. -MAT    *m_transp(in,out)
  2996. -MAT    *in, *out;
  2997. -{
  2998. -    int    i, j;
  2999. -    int    in_situ;
  3000. -    Real    tmp;
  3001. -
  3002. -    if ( in == (MAT *)NULL )
  3003. -        error(E_NULL,"m_transp");
  3004. -    if ( in == out && in->n != in->m )
  3005. -        error(E_INSITU2,"m_transp");
  3006. -    in_situ = ( in == out );
  3007. -    if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m )
  3008. -        out = m_resize(out,in->n,in->m);
  3009. -
  3010. -    if ( ! in_situ )
  3011. -        for ( i = 0; i < in->m; i++ )
  3012. -            for ( j = 0; j < in->n; j++ )
  3013. -                out->me[j][i] = in->me[i][j];
  3014. -    else
  3015. -        for ( i = 1; i < in->m; i++ )
  3016. -            for ( j = 0; j < i; j++ )
  3017. -            {    tmp = in->me[i][j];
  3018. -                in->me[i][j] = in->me[j][i];
  3019. -                in->me[j][i] = tmp;
  3020. -            }
  3021. -
  3022. -    return out;
  3023. -}
  3024. -
  3025. -/* swap_rows -- swaps rows i and j of matrix A upto column lim */
  3026. -MAT    *swap_rows(A,i,j,lo,hi)
  3027. -MAT    *A;
  3028. -int    i, j, lo, hi;
  3029. -{
  3030. -    int    k;
  3031. -    Real    **A_me, tmp;
  3032. -
  3033. -    if ( ! A )
  3034. -        error(E_NULL,"swap_rows");
  3035. -    if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  3036. -        error(E_SIZES,"swap_rows");
  3037. -    lo = max(0,lo);
  3038. -    hi = min(hi,A->n-1);
  3039. -    A_me = A->me;
  3040. -
  3041. -    for ( k = lo; k <= hi; k++ )
  3042. -    {
  3043. -        tmp = A_me[k][i];
  3044. -        A_me[k][i] = A_me[k][j];
  3045. -        A_me[k][j] = tmp;
  3046. -    }
  3047. -    return A;
  3048. -}
  3049. -
  3050. -/* swap_cols -- swap columns i and j of matrix A upto row lim */
  3051. -MAT    *swap_cols(A,i,j,lo,hi)
  3052. -MAT    *A;
  3053. -int    i, j, lo, hi;
  3054. -{
  3055. -    int    k;
  3056. -    Real    **A_me, tmp;
  3057. -
  3058. -    if ( ! A )
  3059. -        error(E_NULL,"swap_cols");
  3060. -    if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  3061. -        error(E_SIZES,"swap_cols");
  3062. -    lo = max(0,lo);
  3063. -    hi = min(hi,A->m-1);
  3064. -    A_me = A->me;
  3065. -
  3066. -    for ( k = lo; k <= hi; k++ )
  3067. -    {
  3068. -        tmp = A_me[i][k];
  3069. -        A_me[i][k] = A_me[j][k];
  3070. -        A_me[j][k] = tmp;
  3071. -    }
  3072. -    return A;
  3073. -}
  3074. -
  3075. -/* ms_mltadd -- matrix-scalar multiply and add
  3076. -    -- may be in situ
  3077. -    -- returns out == A1 + s*A2 */
  3078. -MAT    *ms_mltadd(A1,A2,s,out)
  3079. -MAT    *A1, *A2, *out;
  3080. -double    s;
  3081. -{
  3082. -    /* register Real    *A1_e, *A2_e, *out_e; */
  3083. -    /* register int    j; */
  3084. -    int    i, m, n;
  3085. -
  3086. -    if ( ! A1 || ! A2 )
  3087. -        error(E_NULL,"ms_mltadd");
  3088. -    if ( A1->m != A2->m || A1->n != A2->n )
  3089. -        error(E_SIZES,"ms_mltadd");
  3090. -
  3091. -    if ( s == 0.0 )
  3092. -        return m_copy(A1,out);
  3093. -    if ( s == 1.0 )
  3094. -        return m_add(A1,A2,out);
  3095. -
  3096. -    tracecatch(out = m_copy(A1,out),"ms_mltadd");
  3097. -
  3098. -    m = A1->m;    n = A1->n;
  3099. -    for ( i = 0; i < m; i++ )
  3100. -    {
  3101. -        __mltadd__(out->me[i],A2->me[i],s,(int)n);
  3102. -        /**************************************************
  3103. -        A1_e = A1->me[i];
  3104. -        A2_e = A2->me[i];
  3105. -        out_e = out->me[i];
  3106. -        for ( j = 0; j < n; j++ )
  3107. -            out_e[j] = A1_e[j] + s*A2_e[j];
  3108. -        **************************************************/
  3109. -    }
  3110. -
  3111. -    return out;
  3112. -}
  3113. -
  3114. -/* mv_mltadd -- matrix-vector multiply and add
  3115. -    -- may not be in situ
  3116. -    -- returns out == v1 + alpha*A*v2 */
  3117. -VEC    *mv_mltadd(v1,v2,A,alpha,out)
  3118. -VEC    *v1, *v2, *out;
  3119. -MAT    *A;
  3120. -double    alpha;
  3121. -{
  3122. -    /* register    int    j; */
  3123. -    int    i, m, n;
  3124. -    Real    *v2_ve, *out_ve;
  3125. -
  3126. -    if ( ! v1 || ! v2 || ! A )
  3127. -        error(E_NULL,"mv_mltadd");
  3128. -    if ( out == v2 )
  3129. -        error(E_INSITU,"mv_mltadd");
  3130. -    if ( v1->dim != A->m || v2->dim != A-> n )
  3131. -        error(E_SIZES,"mv_mltadd");
  3132. -
  3133. -    tracecatch(out = v_copy(v1,out),"mv_mltadd");
  3134. -
  3135. -    v2_ve = v2->ve;    out_ve = out->ve;
  3136. -    m = A->m;    n = A->n;
  3137. -
  3138. -    if ( alpha == 0.0 )
  3139. -        return out;
  3140. -
  3141. -    for ( i = 0; i < m; i++ )
  3142. -    {
  3143. -        out_ve[i] += alpha*__ip__(A->me[i],v2_ve,(int)n);
  3144. -        /**************************************************
  3145. -        A_e = A->me[i];
  3146. -        sum = 0.0;
  3147. -        for ( j = 0; j < n; j++ )
  3148. -            sum += A_e[j]*v2_ve[j];
  3149. -        out_ve[i] = v1->ve[i] + alpha*sum;
  3150. -        **************************************************/
  3151. -    }
  3152. -
  3153. -    return out;
  3154. -}
  3155. -
  3156. -/* vm_mltadd -- vector-matrix multiply and add
  3157. -    -- may not be in situ
  3158. -    -- returns out' == v1' + v2'*A */
  3159. -VEC    *vm_mltadd(v1,v2,A,alpha,out)
  3160. -VEC    *v1, *v2, *out;
  3161. -MAT    *A;
  3162. -double    alpha;
  3163. -{
  3164. -    int    /* i, */ j, m, n;
  3165. -    Real    tmp, /* *A_e, */ *out_ve;
  3166. -
  3167. -    if ( ! v1 || ! v2 || ! A )
  3168. -        error(E_NULL,"vm_mltadd");
  3169. -    if ( v2 == out )
  3170. -        error(E_INSITU,"vm_mltadd");
  3171. -    if ( v1->dim != A->n || A->m != v2->dim )
  3172. -        error(E_SIZES,"vm_mltadd");
  3173. -
  3174. -    tracecatch(out = v_copy(v1,out),"vm_mltadd");
  3175. -
  3176. -    out_ve = out->ve;    m = A->m;    n = A->n;
  3177. -    for ( j = 0; j < m; j++ )
  3178. -    {
  3179. -        tmp = v2->ve[j]*alpha;
  3180. -        if ( tmp != 0.0 )
  3181. -            __mltadd__(out_ve,A->me[j],tmp,(int)n);
  3182. -        /**************************************************
  3183. -        A_e = A->me[j];
  3184. -        for ( i = 0; i < n; i++ )
  3185. -            out_ve[i] += A_e[i]*tmp;
  3186. -        **************************************************/
  3187. -    }
  3188. -
  3189. -    return out;
  3190. -}
  3191. -
  3192. //GO.SYSIN DD matop.c
  3193. echo pxop.c 1>&2
  3194. sed >pxop.c <<'//GO.SYSIN DD pxop.c' 's/^-//'
  3195. -
  3196. -/**************************************************************************
  3197. -**
  3198. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3199. -**
  3200. -**                 Meschach Library
  3201. -** 
  3202. -** This Meschach Library is provided "as is" without any express 
  3203. -** or implied warranty of any kind with respect to this software. 
  3204. -** In particular the authors shall not be liable for any direct, 
  3205. -** indirect, special, incidental or consequential damages arising 
  3206. -** in any way from use of the software.
  3207. -** 
  3208. -** Everyone is granted permission to copy, modify and redistribute this
  3209. -** Meschach Library, provided:
  3210. -**  1.  All copies contain this copyright notice.
  3211. -**  2.  All modified copies shall carry a notice stating who
  3212. -**      made the last modification and the date of such modification.
  3213. -**  3.  No charge is made for this software or works derived from it.  
  3214. -**      This clause shall not be construed as constraining other software
  3215. -**      distributed on the same medium as this software, nor is a
  3216. -**      distribution fee considered a charge.
  3217. -**
  3218. -***************************************************************************/
  3219. -
  3220. -
  3221. -/* pxop.c 1.5 12/03/87 */
  3222. -
  3223. -
  3224. -#include    <stdio.h>
  3225. -#include    "matrix.h"
  3226. -
  3227. -static    char    rcsid[] = "$Id: pxop.c,v 1.5 1994/03/23 23:58:50 des Exp $";
  3228. -
  3229. -/**********************************************************************
  3230. -Note: A permutation is often interpreted as a matrix
  3231. -        (i.e. a permutation matrix).
  3232. -    A permutation px represents a permutation matrix P where
  3233. -        P[i][j] == 1 if and only if px->pe[i] == j
  3234. -**********************************************************************/
  3235. -
  3236. -
  3237. -/* px_inv -- invert permutation -- in situ
  3238. -    -- taken from ACM Collected Algorithms #250 */
  3239. -PERM    *px_inv(px,out)
  3240. -PERM    *px, *out;
  3241. -{
  3242. -    int    i, j, k, n, *p;
  3243. -    
  3244. -    out = px_copy(px, out);
  3245. -    n = out->size;
  3246. -    p = (int *)(out->pe);
  3247. -    for ( n--; n>=0; n-- )
  3248. -    {
  3249. -    i = p[n];
  3250. -    if ( i < 0 )    p[n] = -1 - i;
  3251. -    else if ( i != n )
  3252. -    {
  3253. -        k = n;
  3254. -        while (TRUE)
  3255. -        {
  3256. -        if ( i < 0 || i >= out->size )
  3257. -            error(E_BOUNDS,"px_inv");
  3258. -        j = p[i];    p[i] = -1 - k;
  3259. -        if ( j == n )
  3260. -        {    p[n] = i;    break;        }
  3261. -        k = i;        i = j;
  3262. -        }
  3263. -    }
  3264. -    }
  3265. -    return out;
  3266. -}
  3267. -
  3268. -/* px_mlt -- permutation multiplication (composition) */
  3269. -PERM    *px_mlt(px1,px2,out)
  3270. -PERM    *px1,*px2,*out;
  3271. -{
  3272. -    u_int    i,size;
  3273. -    
  3274. -    if ( px1==(PERM *)NULL || px2==(PERM *)NULL )
  3275. -    error(E_NULL,"px_mlt");
  3276. -    if ( px1->size != px2->size )
  3277. -    error(E_SIZES,"px_mlt");
  3278. -    if ( px1 == out || px2 == out )
  3279. -    error(E_INSITU,"px_mlt");
  3280. -    if ( out==(PERM *)NULL || out->size < px1->size )
  3281. -    out = px_resize(out,px1->size);
  3282. -    
  3283. -    size = px1->size;
  3284. -    for ( i=0; i<size; i++ )
  3285. -    if ( px2->pe[i] >= size )
  3286. -        error(E_BOUNDS,"px_mlt");
  3287. -    else
  3288. -        out->pe[i] = px1->pe[px2->pe[i]];
  3289. -    
  3290. -    return out;
  3291. -}
  3292. -
  3293. -/* px_vec -- permute vector */
  3294. -VEC    *px_vec(px,vector,out)
  3295. -PERM    *px;
  3296. -VEC    *vector,*out;
  3297. -{
  3298. -    u_int    old_i, i, size, start;
  3299. -    Real    tmp;
  3300. -    
  3301. -    if ( px==(PERM *)NULL || vector==(VEC *)NULL )
  3302. -    error(E_NULL,"px_vec");
  3303. -    if ( px->size > vector->dim )
  3304. -    error(E_SIZES,"px_vec");
  3305. -    if ( out==(VEC *)NULL || out->dim < vector->dim )
  3306. -    out = v_resize(out,vector->dim);
  3307. -    
  3308. -    size = px->size;
  3309. -    if ( size == 0 )
  3310. -    return v_copy(vector,out);
  3311. -    if ( out != vector )
  3312. -    {
  3313. -    for ( i=0; i<size; i++ )
  3314. -        if ( px->pe[i] >= size )
  3315. -        error(E_BOUNDS,"px_vec");
  3316. -        else
  3317. -        out->ve[i] = vector->ve[px->pe[i]];
  3318. -    }
  3319. -    else
  3320. -    {    /* in situ algorithm */
  3321. -    start = 0;
  3322. -    while ( start < size )
  3323. -    {
  3324. -        old_i = start;
  3325. -        i = px->pe[old_i];
  3326. -        if ( i >= size )
  3327. -        {
  3328. -        start++;
  3329. -        continue;
  3330. -        }
  3331. -        tmp = vector->ve[start];
  3332. -        while ( TRUE )
  3333. -        {
  3334. -        vector->ve[old_i] = vector->ve[i];
  3335. -        px->pe[old_i] = i+size;
  3336. -        old_i = i;
  3337. -        i = px->pe[old_i];
  3338. -        if ( i >= size )
  3339. -            break;
  3340. -        if ( i == start )
  3341. -        {
  3342. -            vector->ve[old_i] = tmp;
  3343. -            px->pe[old_i] = i+size;
  3344. -            break;
  3345. -        }
  3346. -        }
  3347. -        start++;
  3348. -    }
  3349. -
  3350. -    for ( i = 0; i < size; i++ )
  3351. -        if ( px->pe[i] < size )
  3352. -        error(E_BOUNDS,"px_vec");
  3353. -        else
  3354. -        px->pe[i] = px->pe[i]-size;
  3355. -    }
  3356. -    
  3357. -    return out;
  3358. -}
  3359. -
  3360. -/* pxinv_vec -- apply the inverse of px to x, returning the result in out */
  3361. -VEC    *pxinv_vec(px,x,out)
  3362. -PERM    *px;
  3363. -VEC    *x, *out;
  3364. -{
  3365. -    u_int    i, size;
  3366. -    
  3367. -    if ( ! px || ! x )
  3368. -    error(E_NULL,"pxinv_vec");
  3369. -    if ( px->size > x->dim )
  3370. -    error(E_SIZES,"pxinv_vec");
  3371. -    /* if ( x == out )
  3372. -    error(E_INSITU,"pxinv_vec"); */
  3373. -    if ( ! out || out->dim < x->dim )
  3374. -    out = v_resize(out,x->dim);
  3375. -    
  3376. -    size = px->size;
  3377. -    if ( size == 0 )
  3378. -    return v_copy(x,out);
  3379. -    if ( out != x )
  3380. -    {
  3381. -    for ( i=0; i<size; i++ )
  3382. -        if ( px->pe[i] >= size )
  3383. -        error(E_BOUNDS,"pxinv_vec");
  3384. -        else
  3385. -        out->ve[px->pe[i]] = x->ve[i];
  3386. -    }
  3387. -    else
  3388. -    {    /* in situ algorithm --- cheat's way out */
  3389. -    px_inv(px,px);
  3390. -    px_vec(px,x,out);
  3391. -    px_inv(px,px);
  3392. -    }
  3393. -
  3394. -    return out;
  3395. -}
  3396. -
  3397. -
  3398. -
  3399. -/* px_transp -- transpose elements of permutation
  3400. -        -- Really multiplying a permutation by a transposition */
  3401. -PERM    *px_transp(px,i1,i2)
  3402. -PERM    *px;        /* permutation to transpose */
  3403. -u_int    i1,i2;        /* elements to transpose */
  3404. -{
  3405. -    u_int    temp;
  3406. -
  3407. -    if ( px==(PERM *)NULL )
  3408. -        error(E_NULL,"px_transp");
  3409. -
  3410. -    if ( i1 < px->size && i2 < px->size )
  3411. -    {
  3412. -        temp = px->pe[i1];
  3413. -        px->pe[i1] = px->pe[i2];
  3414. -        px->pe[i2] = temp;
  3415. -    }
  3416. -
  3417. -    return px;
  3418. -}
  3419. -
  3420. -/* myqsort -- a cheap implementation of Quicksort on integers
  3421. -        -- returns number of swaps */
  3422. -static int myqsort(a,num)
  3423. -int    *a, num;
  3424. -{
  3425. -    int    i, j, tmp, v;
  3426. -    int    numswaps;
  3427. -
  3428. -    numswaps = 0;
  3429. -    if ( num <= 1 )
  3430. -        return 0;
  3431. -
  3432. -    i = 0;    j = num;    v = a[0];
  3433. -    for ( ; ; )
  3434. -    {
  3435. -        while ( a[++i] < v )
  3436. -            ;
  3437. -        while ( a[--j] > v )
  3438. -            ;
  3439. -        if ( i >= j )    break;
  3440. -
  3441. -        tmp = a[i];
  3442. -        a[i] = a[j];
  3443. -        a[j] = tmp;
  3444. -        numswaps++;
  3445. -    }
  3446. -
  3447. -    tmp = a[0];
  3448. -    a[0] = a[j];
  3449. -    a[j] = tmp;
  3450. -    if ( j != 0 )
  3451. -        numswaps++;
  3452. -
  3453. -    numswaps += myqsort(&a[0],j);
  3454. -    numswaps += myqsort(&a[j+1],num-(j+1));
  3455. -
  3456. -    return numswaps;
  3457. -}
  3458. -
  3459. -
  3460. -/* px_sign -- compute the ``sign'' of a permutation = +/-1 where
  3461. -        px is the product of an even/odd # transpositions */
  3462. -int    px_sign(px)
  3463. -PERM    *px;
  3464. -{
  3465. -    int    numtransp;
  3466. -    PERM    *px2;
  3467. -
  3468. -    if ( px==(PERM *)NULL )
  3469. -        error(E_NULL,"px_sign");
  3470. -    px2 = px_copy(px,PNULL);
  3471. -    numtransp = myqsort(px2->pe,px2->size);
  3472. -    px_free(px2);
  3473. -
  3474. -    return ( numtransp % 2 ) ? -1 : 1;
  3475. -}
  3476. -
  3477. -
  3478. -/* px_cols -- permute columns of matrix A; out = A.px'
  3479. -    -- May NOT be in situ */
  3480. -MAT    *px_cols(px,A,out)
  3481. -PERM    *px;
  3482. -MAT    *A, *out;
  3483. -{
  3484. -    int    i, j, m, n, px_j;
  3485. -    Real    **A_me, **out_me;
  3486. -#ifdef ANSI_C
  3487. -    MAT    *m_get(int, int);
  3488. -#else
  3489. -    extern MAT    *m_get();
  3490. -#endif
  3491. -
  3492. -    if ( ! A || ! px )
  3493. -        error(E_NULL,"px_cols");
  3494. -    if ( px->size != A->n )
  3495. -        error(E_SIZES,"px_cols");
  3496. -    if ( A == out )
  3497. -        error(E_INSITU,"px_cols");
  3498. -    m = A->m;    n = A->n;
  3499. -    if ( ! out || out->m != m || out->n != n )
  3500. -        out = m_get(m,n);
  3501. -    A_me = A->me;    out_me = out->me;
  3502. -
  3503. -    for ( j = 0; j < n; j++ )
  3504. -    {
  3505. -        px_j = px->pe[j];
  3506. -        if ( px_j >= n )
  3507. -            error(E_BOUNDS,"px_cols");
  3508. -        for ( i = 0; i < m; i++ )
  3509. -            out_me[i][px_j] = A_me[i][j];
  3510. -    }
  3511. -
  3512. -    return out;
  3513. -}
  3514. -
  3515. -/* px_rows -- permute columns of matrix A; out = px.A
  3516. -    -- May NOT be in situ */
  3517. -MAT    *px_rows(px,A,out)
  3518. -PERM    *px;
  3519. -MAT    *A, *out;
  3520. -{
  3521. -    int    i, j, m, n, px_i;
  3522. -    Real    **A_me, **out_me;
  3523. -#ifdef ANSI_C
  3524. -    MAT    *m_get(int, int);
  3525. -#else
  3526. -    extern MAT    *m_get();
  3527. -#endif
  3528. -
  3529. -    if ( ! A || ! px )
  3530. -        error(E_NULL,"px_rows");
  3531. -    if ( px->size != A->m )
  3532. -        error(E_SIZES,"px_rows");
  3533. -    if ( A == out )
  3534. -        error(E_INSITU,"px_rows");
  3535. -    m = A->m;    n = A->n;
  3536. -    if ( ! out || out->m != m || out->n != n )
  3537. -        out = m_get(m,n);
  3538. -    A_me = A->me;    out_me = out->me;
  3539. -
  3540. -    for ( i = 0; i < m; i++ )
  3541. -    {
  3542. -        px_i = px->pe[i];
  3543. -        if ( px_i >= m )
  3544. -            error(E_BOUNDS,"px_rows");
  3545. -        for ( j = 0; j < n; j++ )
  3546. -            out_me[i][j] = A_me[px_i][j];
  3547. -    }
  3548. -
  3549. -    return out;
  3550. -}
  3551. -
  3552. //GO.SYSIN DD pxop.c
  3553. echo submat.c 1>&2
  3554. sed >submat.c <<'//GO.SYSIN DD submat.c' 's/^-//'
  3555. -
  3556. -/**************************************************************************
  3557. -**
  3558. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3559. -**
  3560. -**                 Meschach Library
  3561. -** 
  3562. -** This Meschach Library is provided "as is" without any express 
  3563. -** or implied warranty of any kind with respect to this software. 
  3564. -** In particular the authors shall not be liable for any direct, 
  3565. -** indirect, special, incidental or consequential damages arising 
  3566. -** in any way from use of the software.
  3567. -** 
  3568. -** Everyone is granted permission to copy, modify and redistribute this
  3569. -** Meschach Library, provided:
  3570. -**  1.  All copies contain this copyright notice.
  3571. -**  2.  All modified copies shall carry a notice stating who
  3572. -**      made the last modification and the date of such modification.
  3573. -**  3.  No charge is made for this software or works derived from it.  
  3574. -**      This clause shall not be construed as constraining other software
  3575. -**      distributed on the same medium as this software, nor is a
  3576. -**      distribution fee considered a charge.
  3577. -**
  3578. -***************************************************************************/
  3579. -
  3580. -
  3581. -/* 1.2 submat.c 11/25/87 */
  3582. -
  3583. -#include    <stdio.h>
  3584. -#include    "matrix.h"
  3585. -
  3586. -static    char    rcsid[] = "$Id: submat.c,v 1.2 1994/01/13 05:28:12 des Exp $";
  3587. -
  3588. -
  3589. -/* get_col -- gets a specified column of a matrix and retruns it as a vector */
  3590. -VEC    *get_col(mat,col,vec)
  3591. -u_int    col;
  3592. -MAT    *mat;
  3593. -VEC    *vec;
  3594. -{
  3595. -   u_int    i;
  3596. -   
  3597. -   if ( mat==(MAT *)NULL )
  3598. -     error(E_NULL,"get_col");
  3599. -   if ( col >= mat->n )
  3600. -     error(E_RANGE,"get_col");
  3601. -   if ( vec==(VEC *)NULL || vec->dim<mat->m )
  3602. -     vec = v_resize(vec,mat->m);
  3603. -   
  3604. -   for ( i=0; i<mat->m; i++ )
  3605. -     vec->ve[i] = mat->me[i][col];
  3606. -   
  3607. -   return (vec);
  3608. -}
  3609. -
  3610. -/* get_row -- gets a specified row of a matrix and retruns it as a vector */
  3611. -VEC    *get_row(mat,row,vec)
  3612. -u_int    row;
  3613. -MAT    *mat;
  3614. -VEC    *vec;
  3615. -{
  3616. -   u_int    i;
  3617. -   
  3618. -   if ( mat==(MAT *)NULL )
  3619. -     error(E_NULL,"get_row");
  3620. -   if ( row >= mat->m )
  3621. -     error(E_RANGE,"get_row");
  3622. -   if ( vec==(VEC *)NULL || vec->dim<mat->n )
  3623. -     vec = v_resize(vec,mat->n);
  3624. -   
  3625. -   for ( i=0; i<mat->n; i++ )
  3626. -     vec->ve[i] = mat->me[row][i];
  3627. -   
  3628. -   return (vec);
  3629. -}
  3630. -
  3631. -/* _set_col -- sets column of matrix to values given in vec (in situ) */
  3632. -MAT    *_set_col(mat,col,vec,i0)
  3633. -MAT    *mat;
  3634. -VEC    *vec;
  3635. -u_int    col,i0;
  3636. -{
  3637. -   u_int    i,lim;
  3638. -   
  3639. -   if ( mat==(MAT *)NULL || vec==(VEC *)NULL )
  3640. -     error(E_NULL,"_set_col");
  3641. -   if ( col >= mat->n )
  3642. -     error(E_RANGE,"_set_col");
  3643. -   lim = min(mat->m,vec->dim);
  3644. -   for ( i=i0; i<lim; i++ )
  3645. -     mat->me[i][col] = vec->ve[i];
  3646. -   
  3647. -   return (mat);
  3648. -}
  3649. -
  3650. -/* _set_row -- sets row of matrix to values given in vec (in situ) */
  3651. -MAT    *_set_row(mat,row,vec,j0)
  3652. -MAT    *mat;
  3653. -VEC    *vec;
  3654. -u_int    row,j0;
  3655. -{
  3656. -   u_int    j,lim;
  3657. -   
  3658. -   if ( mat==(MAT *)NULL || vec==(VEC *)NULL )
  3659. -     error(E_NULL,"_set_row");
  3660. -   if ( row >= mat->m )
  3661. -     error(E_RANGE,"_set_row");
  3662. -   lim = min(mat->n,vec->dim);
  3663. -   for ( j=j0; j<lim; j++ )
  3664. -     mat->me[row][j] = vec->ve[j];
  3665. -   
  3666. -   return (mat);
  3667. -}
  3668. -
  3669. -/* sub_mat -- returns sub-matrix of old which is formed by the rectangle
  3670. -   from (row1,col1) to (row2,col2)
  3671. -   -- Note: storage is shared so that altering the "new"
  3672. -   matrix will alter the "old" matrix */
  3673. -MAT    *sub_mat(old,row1,col1,row2,col2,new)
  3674. -MAT    *old,*new;
  3675. -u_int    row1,col1,row2,col2;
  3676. -{
  3677. -   u_int    i;
  3678. -   
  3679. -   if ( old==(MAT *)NULL )
  3680. -     error(E_NULL,"sub_mat");
  3681. -   if ( row1 > row2 || col1 > col2 || row2 >= old->m || col2 >= old->n )
  3682. -     error(E_RANGE,"sub_mat");
  3683. -   if ( new==(MAT *)NULL || new->m < row2-row1+1 )
  3684. -   {
  3685. -      new = NEW(MAT);
  3686. -      new->me = NEW_A(row2-row1+1,Real *);
  3687. -      if ( new==(MAT *)NULL || new->me==(Real **)NULL )
  3688. -    error(E_MEM,"sub_mat");
  3689. -      else if (mem_info_is_on()) {
  3690. -     mem_bytes(TYPE_MAT,0,sizeof(MAT)+
  3691. -              (row2-row1+1)*sizeof(Real *));
  3692. -      }
  3693. -      
  3694. -   }
  3695. -   new->m = row2-row1+1;
  3696. -   
  3697. -   new->n = col2-col1+1;
  3698. -   
  3699. -   new->base = (Real *)NULL;
  3700. -   
  3701. -   for ( i=0; i < new->m; i++ )
  3702. -     new->me[i] = (old->me[i+row1]) + col1;
  3703. -   
  3704. -   return (new);
  3705. -}
  3706. -
  3707. -
  3708. -/* sub_vec -- returns sub-vector which is formed by the elements i1 to i2
  3709. -   -- as for sub_mat, storage is shared */
  3710. -VEC    *sub_vec(old,i1,i2,new)
  3711. -VEC    *old, *new;
  3712. -int    i1, i2;
  3713. -{
  3714. -   if ( old == (VEC *)NULL )
  3715. -     error(E_NULL,"sub_vec");
  3716. -   if ( i1 > i2 || old->dim < i2 )
  3717. -     error(E_RANGE,"sub_vec");
  3718. -   
  3719. -   if ( new == (VEC *)NULL )
  3720. -     new = NEW(VEC);
  3721. -   if ( new == (VEC *)NULL )
  3722. -     error(E_MEM,"sub_vec");
  3723. -   else if (mem_info_is_on()) {
  3724. -      mem_bytes(TYPE_VEC,0,sizeof(VEC));
  3725. -   }
  3726. -   
  3727. -   
  3728. -   new->dim = i2 - i1 + 1;
  3729. -   new->ve = &(old->ve[i1]);
  3730. -   
  3731. -   return new;
  3732. -}
  3733. //GO.SYSIN DD submat.c
  3734. echo init.c 1>&2
  3735. sed >init.c <<'//GO.SYSIN DD init.c' 's/^-//'
  3736. -
  3737. -/**************************************************************************
  3738. -**
  3739. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3740. -**
  3741. -**                 Meschach Library
  3742. -** 
  3743. -** This Meschach Library is provided "as is" without any express 
  3744. -** or implied warranty of any kind with respect to this software. 
  3745. -** In particular the authors shall not be liable for any direct, 
  3746. -** indirect, special, incidental or consequential damages arising 
  3747. -** in any way from use of the software.
  3748. -** 
  3749. -** Everyone is granted permission to copy, modify and redistribute this
  3750. -** Meschach Library, provided:
  3751. -**  1.  All copies contain this copyright notice.
  3752. -**  2.  All modified copies shall carry a notice stating who
  3753. -**      made the last modification and the date of such modification.
  3754. -**  3.  No charge is made for this software or works derived from it.  
  3755. -**      This clause shall not be construed as constraining other software
  3756. -**      distributed on the same medium as this software, nor is a
  3757. -**      distribution fee considered a charge.
  3758. -**
  3759. -***************************************************************************/
  3760. -
  3761. -
  3762. -/*
  3763. -    This is a file of routines for zero-ing, and initialising
  3764. -    vectors, matrices and permutations.
  3765. -    This is to be included in the matrix.a library
  3766. -*/
  3767. -
  3768. -static    char    rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";
  3769. -
  3770. -#include    <stdio.h>
  3771. -#include    "matrix.h"
  3772. -
  3773. -/* v_zero -- zero the vector x */
  3774. -VEC    *v_zero(x)
  3775. -VEC    *x;
  3776. -{
  3777. -    if ( x == VNULL )
  3778. -        error(E_NULL,"v_zero");
  3779. -
  3780. -    __zero__(x->ve,x->dim);
  3781. -    /* for ( i = 0; i < x->dim; i++ )
  3782. -        x->ve[i] = 0.0; */
  3783. -
  3784. -    return x;
  3785. -}
  3786. -
  3787. -
  3788. -/* iv_zero -- zero the vector ix */
  3789. -IVEC    *iv_zero(ix)
  3790. -IVEC    *ix;
  3791. -{
  3792. -   int i;
  3793. -   
  3794. -   if ( ix == IVNULL )
  3795. -     error(E_NULL,"iv_zero");
  3796. -   
  3797. -   for ( i = 0; i < ix->dim; i++ )
  3798. -     ix->ive[i] = 0; 
  3799. -   
  3800. -   return ix;
  3801. -}
  3802. -
  3803. -
  3804. -/* m_zero -- zero the matrix A */
  3805. -MAT    *m_zero(A)
  3806. -MAT    *A;
  3807. -{
  3808. -    int    i, A_m, A_n;
  3809. -    Real    **A_me;
  3810. -
  3811. -    if ( A == MNULL )
  3812. -        error(E_NULL,"m_zero");
  3813. -
  3814. -    A_m = A->m;    A_n = A->n;    A_me = A->me;
  3815. -    for ( i = 0; i < A_m; i++ )
  3816. -        __zero__(A_me[i],A_n);
  3817. -        /* for ( j = 0; j < A_n; j++ )
  3818. -            A_me[i][j] = 0.0; */
  3819. -
  3820. -    return A;
  3821. -}
  3822. -
  3823. -/* mat_id -- set A to being closest to identity matrix as possible
  3824. -    -- i.e. A[i][j] == 1 if i == j and 0 otherwise */
  3825. -MAT    *m_ident(A)
  3826. -MAT    *A;
  3827. -{
  3828. -    int    i, size;
  3829. -
  3830. -    if ( A == MNULL )
  3831. -        error(E_NULL,"m_ident");
  3832. -
  3833. -    m_zero(A);
  3834. -    size = min(A->m,A->n);
  3835. -    for ( i = 0; i < size; i++ )
  3836. -        A->me[i][i] = 1.0;
  3837. -
  3838. -    return A;
  3839. -}
  3840. -    
  3841. -/* px_ident -- set px to identity permutation */
  3842. -PERM    *px_ident(px)
  3843. -PERM    *px;
  3844. -{
  3845. -    int    i, px_size;
  3846. -    u_int    *px_pe;
  3847. -
  3848. -    if ( px == PNULL )
  3849. -        error(E_NULL,"px_ident");
  3850. -
  3851. -    px_size = px->size;    px_pe = px->pe;
  3852. -    for ( i = 0; i < px_size; i++ )
  3853. -        px_pe[i] = i;
  3854. -
  3855. -    return px;
  3856. -}
  3857. -
  3858. -/* Pseudo random number generator data structures */
  3859. -/* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
  3860. -   The Art of Computer Programming" sections 3.2-3.3 */
  3861. -
  3862. -#ifdef ANSI_C
  3863. -#ifndef LONG_MAX
  3864. -#include    <limits.h>
  3865. -#endif
  3866. -#endif
  3867. -
  3868. -#ifdef LONG_MAX
  3869. -#define MODULUS    LONG_MAX
  3870. -#else
  3871. -#define MODULUS    1000000000L    /* assuming long's at least 32 bits long */
  3872. -#endif
  3873. -#define MZ    0L
  3874. -
  3875. -static long mrand_list[56];
  3876. -static int  started = FALSE;
  3877. -static int  inext = 0, inextp = 31;
  3878. -
  3879. -
  3880. -/* mrand -- pseudo-random number generator */
  3881. -#ifdef ANSI_C
  3882. -double mrand(void)
  3883. -#else
  3884. -double mrand()
  3885. -#endif
  3886. -{
  3887. -    long    lval;
  3888. -    static Real  factor = 1.0/((Real)MODULUS);
  3889. -    
  3890. -    if ( ! started )
  3891. -    smrand(3127);
  3892. -    
  3893. -    inext = (inext >= 54) ? 0 : inext+1;
  3894. -    inextp = (inextp >= 54) ? 0 : inextp+1;
  3895. -
  3896. -    lval = mrand_list[inext]-mrand_list[inextp];
  3897. -    if ( lval < 0L )
  3898. -    lval += MODULUS;
  3899. -    mrand_list[inext] = lval;
  3900. -    
  3901. -    return (double)lval*factor;
  3902. -}
  3903. -
  3904. -/* mrandlist -- fills the array a[] with len random numbers */
  3905. -void    mrandlist(a, len)
  3906. -Real    a[];
  3907. -int    len;
  3908. -{
  3909. -    int        i;
  3910. -    long    lval;
  3911. -    static Real  factor = 1.0/((Real)MODULUS);
  3912. -    
  3913. -    if ( ! started )
  3914. -    smrand(3127);
  3915. -    
  3916. -    for ( i = 0; i < len; i++ )
  3917. -    {
  3918. -    inext = (inext >= 54) ? 0 : inext+1;
  3919. -    inextp = (inextp >= 54) ? 0 : inextp+1;
  3920. -    
  3921. -    lval = mrand_list[inext]-mrand_list[inextp];
  3922. -    if ( lval < 0L )
  3923. -        lval += MODULUS;
  3924. -    mrand_list[inext] = lval;
  3925. -    
  3926. -    a[i] = (Real)lval*factor;
  3927. -    }
  3928. -}
  3929. -
  3930. -
  3931. -/* smrand -- set seed for mrand() */
  3932. -void smrand(seed)
  3933. -int    seed;
  3934. -{
  3935. -    int        i;
  3936. -
  3937. -    mrand_list[0] = (123413*seed) % MODULUS;
  3938. -    for ( i = 1; i < 55; i++ )
  3939. -    mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;
  3940. -
  3941. -    started = TRUE;
  3942. -
  3943. -    /* run mrand() through the list sufficient times to
  3944. -       thoroughly randomise the array */
  3945. -    for ( i = 0; i < 55*55; i++ )
  3946. -    mrand();
  3947. -}
  3948. -#undef MODULUS
  3949. -#undef MZ
  3950. -#undef FAC
  3951. -
  3952. -/* v_rand -- initialises x to be a random vector, components
  3953. -    independently & uniformly ditributed between 0 and 1 */
  3954. -VEC    *v_rand(x)
  3955. -VEC    *x;
  3956. -{
  3957. -    /* int    i; */
  3958. -
  3959. -    if ( ! x )
  3960. -        error(E_NULL,"v_rand");
  3961. -
  3962. -    /* for ( i = 0; i < x->dim; i++ ) */
  3963. -        /* x->ve[i] = rand()/((Real)MAX_RAND); */
  3964. -        /* x->ve[i] = mrand(); */
  3965. -    mrandlist(x->ve,x->dim);
  3966. -
  3967. -    return x;
  3968. -}
  3969. -
  3970. -/* m_rand -- initialises A to be a random vector, components
  3971. -    independently & uniformly distributed between 0 and 1 */
  3972. -MAT    *m_rand(A)
  3973. -MAT    *A;
  3974. -{
  3975. -    int    i /* , j */;
  3976. -
  3977. -    if ( ! A )
  3978. -        error(E_NULL,"m_rand");
  3979. -
  3980. -    for ( i = 0; i < A->m; i++ )
  3981. -        /* for ( j = 0; j < A->n; j++ ) */
  3982. -            /* A->me[i][j] = rand()/((Real)MAX_RAND); */
  3983. -            /* A->me[i][j] = mrand(); */
  3984. -        mrandlist(A->me[i],A->n);
  3985. -
  3986. -    return A;
  3987. -}
  3988. -
  3989. -/* v_ones -- fills x with one's */
  3990. -VEC    *v_ones(x)
  3991. -VEC    *x;
  3992. -{
  3993. -    int    i;
  3994. -
  3995. -    if ( ! x )
  3996. -        error(E_NULL,"v_ones");
  3997. -
  3998. -    for ( i = 0; i < x->dim; i++ )
  3999. -        x->ve[i] = 1.0;
  4000. -
  4001. -    return x;
  4002. -}
  4003. -
  4004. -/* m_ones -- fills matrix with one's */
  4005. -MAT    *m_ones(A)
  4006. -MAT    *A;
  4007. -{
  4008. -    int    i, j;
  4009. -
  4010. -    if ( ! A )
  4011. -        error(E_NULL,"m_ones");
  4012. -
  4013. -    for ( i = 0; i < A->m; i++ )
  4014. -        for ( j = 0; j < A->n; j++ )
  4015. -            A->me[i][j] = 1.0;
  4016. -
  4017. -    return A;
  4018. -}
  4019. -
  4020. -/* v_count -- initialises x so that x->ve[i] == i */
  4021. -VEC    *v_count(x)
  4022. -VEC    *x;
  4023. -{
  4024. -    int    i;
  4025. -
  4026. -    if ( ! x )
  4027. -        error(E_NULL,"v_count");
  4028. -
  4029. -    for ( i = 0; i < x->dim; i++ )
  4030. -        x->ve[i] = (Real)i;
  4031. -
  4032. -    return x;
  4033. -}
  4034. //GO.SYSIN DD init.c
  4035. echo otherio.c 1>&2
  4036. sed >otherio.c <<'//GO.SYSIN DD otherio.c' 's/^-//'
  4037. -
  4038. -/**************************************************************************
  4039. -**
  4040. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4041. -**
  4042. -**                 Meschach Library
  4043. -** 
  4044. -** This Meschach Library is provided "as is" without any express 
  4045. -** or implied warranty of any kind with respect to this software. 
  4046. -** In particular the authors shall not be liable for any direct, 
  4047. -** indirect, special, incidental or consequential damages arising 
  4048. -** in any way from use of the software.
  4049. -** 
  4050. -** Everyone is granted permission to copy, modify and redistribute this
  4051. -** Meschach Library, provided:
  4052. -**  1.  All copies contain this copyright notice.
  4053. -**  2.  All modified copies shall carry a notice stating who
  4054. -**      made the last modification and the date of such modification.
  4055. -**  3.  No charge is made for this software or works derived from it.  
  4056. -**      This clause shall not be construed as constraining other software
  4057. -**      distributed on the same medium as this software, nor is a
  4058. -**      distribution fee considered a charge.
  4059. -**
  4060. -***************************************************************************/
  4061. -
  4062. -
  4063. -/*
  4064. -    File for doing assorted I/O operations not invlolving
  4065. -    MAT/VEC/PERM objects
  4066. -*/
  4067. -static    char    rcsid[] = "$Id: otherio.c,v 1.2 1994/01/13 05:34:52 des Exp $";
  4068. -
  4069. -#include    <stdio.h>
  4070. -#include    <ctype.h>
  4071. -#include    "matrix.h"
  4072. -
  4073. -
  4074. -
  4075. -/* scratch area -- enough for a single line */
  4076. -static    char    scratch[MAXLINE+1];
  4077. -
  4078. -/* default value for fy_or_n */
  4079. -static    int    y_n_dflt = TRUE;
  4080. -
  4081. -/* fy_or_n -- yes-or-no to question is string s
  4082. -    -- question written to stderr, input from fp 
  4083. -    -- if fp is NOT a tty then return y_n_dflt */
  4084. -int    fy_or_n(fp,s)
  4085. -FILE    *fp;
  4086. -char    *s;
  4087. -{
  4088. -    char    *cp;
  4089. -
  4090. -    if ( ! isatty(fileno(fp)) )
  4091. -        return y_n_dflt;
  4092. -
  4093. -    for ( ; ; )
  4094. -    {
  4095. -        fprintf(stderr,"%s (y/n) ? ",s);
  4096. -        if ( fgets(scratch,MAXLINE,fp)==NULL )
  4097. -            error(E_INPUT,"fy_or_n");
  4098. -        cp = scratch;
  4099. -        while ( isspace(*cp) )
  4100. -            cp++;
  4101. -        if ( *cp == 'y' || *cp == 'Y' )
  4102. -            return TRUE;
  4103. -        if ( *cp == 'n' || *cp == 'N' )
  4104. -            return FALSE;
  4105. -        fprintf(stderr,"Please reply with 'y' or 'Y' for yes ");
  4106. -        fprintf(stderr,"and 'n' or 'N' for no.\n");
  4107. -    }
  4108. -}
  4109. -
  4110. -/* yn_dflt -- sets the value of y_n_dflt to val */
  4111. -int    yn_dflt(val)
  4112. -int    val;
  4113. -{    return y_n_dflt = val;        }
  4114. -
  4115. -/* fin_int -- return integer read from file/stream fp
  4116. -    -- prompt s on stderr if fp is a tty
  4117. -    -- check that x lies between low and high: re-prompt if
  4118. -        fp is a tty, error exit otherwise
  4119. -    -- ignore check if low > high        */
  4120. -int    fin_int(fp,s,low,high)
  4121. -FILE    *fp;
  4122. -char    *s;
  4123. -int    low, high;
  4124. -{
  4125. -    int    retcode, x;
  4126. -
  4127. -    if ( ! isatty(fileno(fp)) )
  4128. -    {
  4129. -        skipjunk(fp);
  4130. -        if ( (retcode=fscanf(fp,"%d",&x)) == EOF )
  4131. -            error(E_INPUT,"fin_int");
  4132. -        if ( retcode <= 0 )
  4133. -            error(E_FORMAT,"fin_int");
  4134. -        if ( low <= high && ( x < low || x > high ) )
  4135. -            error(E_BOUNDS,"fin_int");
  4136. -        return x;
  4137. -    }
  4138. -
  4139. -    for ( ; ; )
  4140. -    {
  4141. -        fprintf(stderr,"%s: ",s);
  4142. -        if ( fgets(scratch,MAXLINE,stdin)==NULL )
  4143. -            error(E_INPUT,"fin_int");
  4144. -        retcode = sscanf(scratch,"%d",&x);
  4145. -        if ( ( retcode==1 && low > high ) ||
  4146. -                    ( x >= low && x <= high ) )
  4147. -            return x;
  4148. -        fprintf(stderr,"Please type an integer in range [%d,%d].\n",
  4149. -                            low,high);
  4150. -    }
  4151. -}
  4152. -
  4153. -
  4154. -/* fin_double -- return double read from file/stream fp
  4155. -    -- prompt s on stderr if fp is a tty
  4156. -    -- check that x lies between low and high: re-prompt if
  4157. -        fp is a tty, error exit otherwise
  4158. -    -- ignore check if low > high        */
  4159. -double    fin_double(fp,s,low,high)
  4160. -FILE    *fp;
  4161. -char    *s;
  4162. -double    low, high;
  4163. -{
  4164. -    Real    retcode, x;
  4165. -
  4166. -    if ( ! isatty(fileno(fp)) )
  4167. -    {
  4168. -        skipjunk(fp);
  4169. -#if REAL == DOUBLE
  4170. -        if ( (retcode=fscanf(fp,"%lf",&x)) == EOF )
  4171. -#elif REAL == FLOAT
  4172. -        if ( (retcode=fscanf(fp,"%f",&x)) == EOF )
  4173. -#endif
  4174. -            error(E_INPUT,"fin_double");
  4175. -        if ( retcode <= 0 )
  4176. -            error(E_FORMAT,"fin_double");
  4177. -        if ( low <= high && ( x < low || x > high ) )
  4178. -            error(E_BOUNDS,"fin_double");
  4179. -        return (double)x;
  4180. -    }
  4181. -
  4182. -    for ( ; ; )
  4183. -    {
  4184. -        fprintf(stderr,"%s: ",s);
  4185. -        if ( fgets(scratch,MAXLINE,stdin)==NULL )
  4186. -            error(E_INPUT,"fin_double");
  4187. -#if REAL == DOUBLE
  4188. -        retcode = sscanf(scratch,"%lf",&x);
  4189. -#elif REAL == FLOAT 
  4190. -        retcode = sscanf(scratch,"%f",&x);
  4191. -#endif
  4192. -        if ( ( retcode==1 && low > high ) ||
  4193. -                    ( x >= low && x <= high ) )
  4194. -            return (double)x;
  4195. -        fprintf(stderr,"Please type an double in range [%g,%g].\n",
  4196. -                            low,high);
  4197. -    }
  4198. -}
  4199. -
  4200. -
  4201. //GO.SYSIN DD otherio.c
  4202. echo machine.c 1>&2
  4203. sed >machine.c <<'//GO.SYSIN DD machine.c' 's/^-//'
  4204. -
  4205. -/**************************************************************************
  4206. -**
  4207. -** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
  4208. -**
  4209. -**                 Meschach Library
  4210. -** 
  4211. -** This Meschach Library is provided "as is" without any express 
  4212. -** or implied warranty of any kind with respect to this software. 
  4213. -** In particular the authors shall not be liable for any direct, 
  4214. -** indirect, special, incidental or consequential damages arising 
  4215. -** in any way from use of the software.
  4216. -** 
  4217. -** Everyone is granted permission to copy, modify and redistribute this
  4218. -** Meschach Library, provided:
  4219. -**  1.  All copies contain this copyright notice.
  4220. -**  2.  All modified copies shall carry a notice stating who
  4221. -**      made the last modification and the date of such modification.
  4222. -**  3.  No charge is made for this software or works derived from it.  
  4223. -**      This clause shall not be construed as constraining other software
  4224. -**      distributed on the same medium as this software, nor is a
  4225. -**      distribution fee considered a charge.
  4226. -**
  4227. -***************************************************************************/
  4228. -
  4229. -/*
  4230. -  This file contains basic routines which are used by the functions
  4231. -  in meschach.a etc.
  4232. -  These are the routines that should be modified in order to take
  4233. -  full advantage of specialised architectures (pipelining, vector
  4234. -  processors etc).
  4235. -  */
  4236. -
  4237. -static    char    *rcsid = "$Id: machine.c,v 1.4 1994/01/13 05:28:56 des Exp $";
  4238. -
  4239. -#include    "machine.h"
  4240. -
  4241. -/* __ip__ -- inner product */
  4242. -double    __ip__(dp1,dp2,len)
  4243. -register Real    *dp1, *dp2;
  4244. -int    len;
  4245. -{
  4246. -#ifdef VUNROLL
  4247. -    register int    len4;
  4248. -    register Real    sum1, sum2, sum3;
  4249. -#endif
  4250. -    register int    i;
  4251. -    register Real     sum;
  4252. -
  4253. -    sum = 0.0;
  4254. -#ifdef VUNROLL
  4255. -    sum1 = sum2 = sum3 = 0.0;
  4256. -    
  4257. -    len4 = len / 4;
  4258. -    len  = len % 4;
  4259. -    
  4260. -    for ( i = 0; i < len4; i++ )
  4261. -    {
  4262. -    sum  += dp1[4*i]*dp2[4*i];
  4263. -    sum1 += dp1[4*i+1]*dp2[4*i+1];
  4264. -    sum2 += dp1[4*i+2]*dp2[4*i+2];
  4265. -    sum3 += dp1[4*i+3]*dp2[4*i+3];
  4266. -    }
  4267. -    sum  += sum1 + sum2 + sum3;
  4268. -    dp1 += 4*len4;    dp2 += 4*len4;
  4269. -#endif
  4270. -    
  4271. -    for ( i = 0; i < len; i++ )
  4272. -    sum  += dp1[i]*dp2[i];
  4273. -    
  4274. -    return sum;
  4275. -}
  4276. -
  4277. -/* __mltadd__ -- scalar multiply and add c.f. v_mltadd() */
  4278. -void    __mltadd__(dp1,dp2,s,len)
  4279. -register Real    *dp1, *dp2;
  4280. -register double s;
  4281. -register int    len;
  4282. -{
  4283. -    register int    i;
  4284. -#ifdef VUNROLL
  4285. -    register int        len4;
  4286. -    
  4287. -    len4 = len / 4;
  4288. -    len  = len % 4;
  4289. -    for ( i = 0; i < len4; i++ )
  4290. -    {
  4291. -    dp1[4*i]   += s*dp2[4*i];
  4292. -    dp1[4*i+1] += s*dp2[4*i+1];
  4293. -    dp1[4*i+2] += s*dp2[4*i+2];
  4294. -    dp1[4*i+3] += s*dp2[4*i+3];
  4295. -    }
  4296. -    dp1 += 4*len4;    dp2 += 4*len4;
  4297. -#endif
  4298. -    
  4299. -    for ( i = 0; i < len; i++ )
  4300. -    dp1[i] += s*dp2[i];
  4301. -}
  4302. -
  4303. -/* __smlt__ scalar multiply array c.f. sv_mlt() */
  4304. -void    __smlt__(dp,s,out,len)
  4305. -register Real    *dp, *out;
  4306. -register double s;
  4307. -register int    len;
  4308. -{
  4309. -    register int    i;
  4310. -    for ( i = 0; i < len; i++ )
  4311. -    out[i] = s*dp[i];
  4312. -}
  4313. -
  4314. -/* __add__ -- add arrays c.f. v_add() */
  4315. -void    __add__(dp1,dp2,out,len)
  4316. -register Real    *dp1, *dp2, *out;
  4317. -register int    len;
  4318. -{
  4319. -    register int    i;
  4320. -    for ( i = 0; i < len; i++ )
  4321. -    out[i] = dp1[i] + dp2[i];
  4322. -}
  4323. -
  4324. -/* __sub__ -- subtract arrays c.f. v_sub() */
  4325. -void    __sub__(dp1,dp2,out,len)
  4326. -register Real    *dp1, *dp2, *out;
  4327. -register int    len;
  4328. -{
  4329. -    register int    i;
  4330. -    for ( i = 0; i < len; i++ )
  4331. -    out[i] = dp1[i] - dp2[i];
  4332. -}
  4333. -
  4334. -/* __zero__ -- zeros an array of floating point numbers */
  4335. -void    __zero__(dp,len)
  4336. -register Real    *dp;
  4337. -register int    len;
  4338. -{
  4339. -#ifdef CHAR0ISDBL0
  4340. -    /* if a floating point zero is equivalent to a string of nulls */
  4341. -    MEM_ZERO((char *)dp,len*sizeof(Real));
  4342. -#else
  4343. -    /* else, need to zero the array entry by entry */
  4344. -    int    i;
  4345. -    for ( i = 0; i < len; i++ )
  4346. -    dp[i] = 0.0;
  4347. -#endif
  4348. -}
  4349. -
  4350. //GO.SYSIN DD machine.c
  4351. echo matlab.c 1>&2
  4352. sed >matlab.c <<'//GO.SYSIN DD matlab.c' 's/^-//'
  4353. -
  4354. -/**************************************************************************
  4355. -**
  4356. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4357. -**
  4358. -**                 Meschach Library
  4359. -** 
  4360. -** This Meschach Library is provided "as is" without any express 
  4361. -** or implied warranty of any kind with respect to this software. 
  4362. -** In particular the authors shall not be liable for any direct, 
  4363. -** indirect, special, incidental or consequential damages arising 
  4364. -** in any way from use of the software.
  4365. -** 
  4366. -** Everyone is granted permission to copy, modify and redistribute this
  4367. -** Meschach Library, provided:
  4368. -**  1.  All copies contain this copyright notice.
  4369. -**  2.  All modified copies shall carry a notice stating who
  4370. -**      made the last modification and the date of such modification.
  4371. -**  3.  No charge is made for this software or works derived from it.  
  4372. -**      This clause shall not be construed as constraining other software
  4373. -**      distributed on the same medium as this software, nor is a
  4374. -**      distribution fee considered a charge.
  4375. -**
  4376. -***************************************************************************/
  4377. -
  4378. -
  4379. -/*
  4380. -    This file contains routines for import/exporting data to/from
  4381. -        MATLAB. The main routines are:
  4382. -            MAT *m_save(FILE *fp,MAT *A,char *name)
  4383. -            VEC *v_save(FILE *fp,VEC *x,char *name)
  4384. -            MAT *m_load(FILE *fp,char **name)
  4385. -*/
  4386. -
  4387. -#include        <stdio.h>
  4388. -#include        "matrix.h"
  4389. -#include    "matlab.h"
  4390. -
  4391. -static char rcsid[] = "$Id: matlab.c,v 1.7 1994/01/13 05:30:09 des Exp $";
  4392. -
  4393. -/* m_save -- save matrix in ".mat" file for MATLAB
  4394. -    -- returns matrix to be saved */
  4395. -MAT     *m_save(fp,A,name)
  4396. -FILE    *fp;
  4397. -MAT     *A;
  4398. -char    *name;
  4399. -{
  4400. -    int     i;
  4401. -    matlab  mat;
  4402. -
  4403. -    if ( ! A )
  4404. -        error(E_NULL,"m_save");
  4405. -
  4406. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  4407. -    mat.m = A->m;
  4408. -    mat.n = A->n;
  4409. -    mat.imag = FALSE;
  4410. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  4411. -
  4412. -    /* write header */
  4413. -    fwrite(&mat,sizeof(matlab),1,fp);
  4414. -    /* write name */
  4415. -    if ( name == (char *)NULL )
  4416. -        fwrite("",sizeof(char),1,fp);
  4417. -    else
  4418. -        fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  4419. -    /* write actual data */
  4420. -    for ( i = 0; i < A->m; i++ )
  4421. -        fwrite(A->me[i],sizeof(Real),(int)(A->n),fp);
  4422. -
  4423. -    return A;
  4424. -}
  4425. -
  4426. -
  4427. -/* v_save -- save vector in ".mat" file for MATLAB
  4428. -    -- saves it as a row vector
  4429. -    -- returns vector to be saved */
  4430. -VEC     *v_save(fp,x,name)
  4431. -FILE    *fp;
  4432. -VEC     *x;
  4433. -char    *name;
  4434. -{
  4435. -    matlab  mat;
  4436. -
  4437. -    if ( ! x )
  4438. -        error(E_NULL,"v_save");
  4439. -
  4440. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  4441. -    mat.m = x->dim;
  4442. -    mat.n = 1;
  4443. -    mat.imag = FALSE;
  4444. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  4445. -
  4446. -    /* write header */
  4447. -    fwrite(&mat,sizeof(matlab),1,fp);
  4448. -    /* write name */
  4449. -    if ( name == (char *)NULL )
  4450. -        fwrite("",sizeof(char),1,fp);
  4451. -    else
  4452. -        fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  4453. -    /* write actual data */
  4454. -    fwrite(x->ve,sizeof(Real),(int)(x->dim),fp);
  4455. -
  4456. -    return x;
  4457. -}
  4458. -
  4459. -/* d_save -- save double in ".mat" file for MATLAB
  4460. -    -- saves it as a row vector
  4461. -    -- returns vector to be saved */
  4462. -double    d_save(fp,x,name)
  4463. -FILE    *fp;
  4464. -double    x;
  4465. -char    *name;
  4466. -{
  4467. -    matlab  mat;
  4468. -    Real x1 = x;
  4469. -
  4470. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  4471. -    mat.m = 1;
  4472. -    mat.n = 1;
  4473. -    mat.imag = FALSE;
  4474. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  4475. -
  4476. -    /* write header */
  4477. -    fwrite(&mat,sizeof(matlab),1,fp);
  4478. -    /* write name */
  4479. -    if ( name == (char *)NULL )
  4480. -        fwrite("",sizeof(char),1,fp);
  4481. -    else
  4482. -        fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  4483. -    /* write actual data */
  4484. -    fwrite(&x1,sizeof(Real),1,fp);
  4485. -
  4486. -    return x;
  4487. -}
  4488. -
  4489. -/* m_load -- loads in a ".mat" file variable as produced by MATLAB
  4490. -    -- matrix returned; imaginary parts ignored */
  4491. -MAT     *m_load(fp,name)
  4492. -FILE    *fp;
  4493. -char    **name;
  4494. -{
  4495. -    MAT     *A;
  4496. -    int     i;
  4497. -    int     m_flag, o_flag, p_flag, t_flag;
  4498. -    float   f_temp;
  4499. -    Real    d_temp;
  4500. -    matlab  mat;
  4501. -
  4502. -    if ( fread(&mat,sizeof(matlab),1,fp) != 1 )
  4503. -        error(E_FORMAT,"m_load");
  4504. -    if ( mat.type >= 10000 )    /* don't load a sparse matrix! */
  4505. -        error(E_FORMAT,"m_load");
  4506. -    m_flag = (mat.type/1000) % 10;
  4507. -    o_flag = (mat.type/100) % 10;
  4508. -    p_flag = (mat.type/10) % 10;
  4509. -    t_flag = (mat.type) % 10;
  4510. -    if ( m_flag != MACH_ID )
  4511. -        error(E_FORMAT,"m_load");
  4512. -    if ( t_flag != 0 )
  4513. -        error(E_FORMAT,"m_load");
  4514. -    if ( p_flag != DOUBLE_PREC && p_flag != SINGLE_PREC )
  4515. -        error(E_FORMAT,"m_load");
  4516. -    *name = (char *)malloc((unsigned)(mat.namlen)+1);
  4517. -    if ( fread(*name,sizeof(char),(unsigned)(mat.namlen),fp) == 0 )
  4518. -        error(E_FORMAT,"m_load");
  4519. -    A = m_get((unsigned)(mat.m),(unsigned)(mat.n));
  4520. -    for ( i = 0; i < A->m*A->n; i++ )
  4521. -    {
  4522. -        if ( p_flag == DOUBLE_PREC )
  4523. -            fread(&d_temp,sizeof(double),1,fp);
  4524. -        else
  4525. -        {
  4526. -            fread(&f_temp,sizeof(float),1,fp);
  4527. -            d_temp = f_temp;
  4528. -        }
  4529. -        if ( o_flag == ROW_ORDER )
  4530. -            A->me[i / A->n][i % A->n] = d_temp;
  4531. -        else if ( o_flag == COL_ORDER )
  4532. -            A->me[i % A->m][i / A->m] = d_temp;
  4533. -        else
  4534. -            error(E_FORMAT,"m_load");
  4535. -    }
  4536. -
  4537. -    if ( mat.imag )         /* skip imaginary part */
  4538. -    for ( i = 0; i < A->m*A->n; i++ )
  4539. -    {
  4540. -        if ( p_flag == DOUBLE_PREC )
  4541. -            fread(&d_temp,sizeof(double),1,fp);
  4542. -        else
  4543. -            fread(&f_temp,sizeof(float),1,fp);
  4544. -    }
  4545. -
  4546. -    return A;
  4547. -}
  4548. -
  4549. //GO.SYSIN DD matlab.c
  4550. echo ivecop.c 1>&2
  4551. sed >ivecop.c <<'//GO.SYSIN DD ivecop.c' 's/^-//'
  4552. -
  4553. -/**************************************************************************
  4554. -**
  4555. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4556. -**
  4557. -**                 Meschach Library
  4558. -** 
  4559. -** This Meschach Library is provided "as is" without any express 
  4560. -** or implied warranty of any kind with respect to this software. 
  4561. -** In particular the authors shall not be liable for any direct, 
  4562. -** indirect, special, incidental or consequential damages arising 
  4563. -** in any way from use of the software.
  4564. -** 
  4565. -** Everyone is granted permission to copy, modify and redistribute this
  4566. -** Meschach Library, provided:
  4567. -**  1.  All copies contain this copyright notice.
  4568. -**  2.  All modified copies shall carry a notice stating who
  4569. -**      made the last modification and the date of such modification.
  4570. -**  3.  No charge is made for this software or works derived from it.  
  4571. -**      This clause shall not be construed as constraining other software
  4572. -**      distributed on the same medium as this software, nor is a
  4573. -**      distribution fee considered a charge.
  4574. -**
  4575. -***************************************************************************/
  4576. -
  4577. -
  4578. -/* ivecop.c  */
  4579. -
  4580. -#include    <stdio.h>
  4581. -#include     "matrix.h"
  4582. -
  4583. -static    char    rcsid[] = "$Id: ivecop.c,v 1.5 1994/01/13 05:45:30 des Exp $";
  4584. -
  4585. -static char    line[MAXLINE];
  4586. -
  4587. -
  4588. -
  4589. -/* iv_get -- get integer vector -- see also memory.c */
  4590. -IVEC    *iv_get(dim)
  4591. -int    dim;
  4592. -{
  4593. -   IVEC    *iv;
  4594. -   /* u_int    i; */
  4595. -   
  4596. -   if (dim < 0)
  4597. -     error(E_NEG,"iv_get");
  4598. -
  4599. -   if ((iv=NEW(IVEC)) == IVNULL )
  4600. -     error(E_MEM,"iv_get");
  4601. -   else if (mem_info_is_on()) {
  4602. -      mem_bytes(TYPE_IVEC,0,sizeof(IVEC));
  4603. -      mem_numvar(TYPE_IVEC,1);
  4604. -   }
  4605. -   
  4606. -   iv->dim = iv->max_dim = dim;
  4607. -   if ((iv->ive = NEW_A(dim,int)) == (int *)NULL )
  4608. -     error(E_MEM,"iv_get");
  4609. -   else if (mem_info_is_on()) {
  4610. -      mem_bytes(TYPE_IVEC,0,dim*sizeof(int));
  4611. -   }
  4612. -   
  4613. -   return (iv);
  4614. -}
  4615. -
  4616. -/* iv_free -- returns iv & asoociated memory back to memory heap */
  4617. -int    iv_free(iv)
  4618. -IVEC    *iv;
  4619. -{
  4620. -   if ( iv==IVNULL || iv->dim > MAXDIM )
  4621. -     /* don't trust it */
  4622. -     return (-1);
  4623. -   
  4624. -   if ( iv->ive == (int *)NULL ) {
  4625. -      if (mem_info_is_on()) {
  4626. -     mem_bytes(TYPE_IVEC,sizeof(IVEC),0);
  4627. -     mem_numvar(TYPE_IVEC,-1);
  4628. -      }
  4629. -      free((char *)iv);
  4630. -   }
  4631. -   else
  4632. -   {
  4633. -      if (mem_info_is_on()) {
  4634. -     mem_bytes(TYPE_IVEC,sizeof(IVEC)+iv->max_dim*sizeof(int),0);
  4635. -     mem_numvar(TYPE_IVEC,-1);
  4636. -      }    
  4637. -      free((char *)iv->ive);
  4638. -      free((char *)iv);
  4639. -   }
  4640. -   
  4641. -   return (0);
  4642. -}
  4643. -
  4644. -/* iv_resize -- returns the IVEC with dimension new_dim
  4645. -   -- iv is set to the zero vector */
  4646. -IVEC    *iv_resize(iv,new_dim)
  4647. -IVEC    *iv;
  4648. -int    new_dim;
  4649. -{
  4650. -   int    i;
  4651. -   
  4652. -   if (new_dim < 0)
  4653. -     error(E_NEG,"iv_resize");
  4654. -
  4655. -   if ( ! iv )
  4656. -     return iv_get(new_dim);
  4657. -   
  4658. -   if (new_dim == iv->dim)
  4659. -     return iv;
  4660. -
  4661. -   if ( new_dim > iv->max_dim )
  4662. -   {
  4663. -      if (mem_info_is_on()) {
  4664. -     mem_bytes(TYPE_IVEC,iv->max_dim*sizeof(int),
  4665. -              new_dim*sizeof(int));
  4666. -      }
  4667. -      iv->ive = RENEW(iv->ive,new_dim,int);
  4668. -      if ( ! iv->ive )
  4669. -    error(E_MEM,"iv_resize");
  4670. -      iv->max_dim = new_dim;
  4671. -   }
  4672. -   if ( iv->dim <= new_dim )
  4673. -     for ( i = iv->dim; i < new_dim; i++ )
  4674. -       iv->ive[i] = 0;
  4675. -   iv->dim = new_dim;
  4676. -   
  4677. -   return iv;
  4678. -}
  4679. -
  4680. -/* iv_copy -- copy integer vector in to out
  4681. -   -- out created/resized if necessary */
  4682. -IVEC    *iv_copy(in,out)
  4683. -IVEC    *in, *out;
  4684. -{
  4685. -   int        i;
  4686. -   
  4687. -   if ( ! in )
  4688. -     error(E_NULL,"iv_copy");
  4689. -   out = iv_resize(out,in->dim);
  4690. -   for ( i = 0; i < in->dim; i++ )
  4691. -     out->ive[i] = in->ive[i];
  4692. -   
  4693. -   return out;
  4694. -}
  4695. -
  4696. -/* iv_move -- move selected pieces of an IVEC
  4697. -    -- moves the length dim0 subvector with initial index i0
  4698. -       to the corresponding subvector of out with initial index i1
  4699. -    -- out is resized if necessary */
  4700. -IVEC    *iv_move(in,i0,dim0,out,i1)
  4701. -IVEC    *in, *out;
  4702. -int    i0, dim0, i1;
  4703. -{
  4704. -    if ( ! in )
  4705. -    error(E_NULL,"iv_move");
  4706. -    if ( i0 < 0 || dim0 < 0 || i1 < 0 ||
  4707. -     i0+dim0 > in->dim )
  4708. -    error(E_BOUNDS,"iv_move");
  4709. -
  4710. -    if ( (! out) || i1+dim0 > out->dim )
  4711. -    out = iv_resize(out,i1+dim0);
  4712. -
  4713. -    MEM_COPY(&(in->ive[i0]),&(out->ive[i1]),dim0*sizeof(int));
  4714. -
  4715. -    return out;
  4716. -}
  4717. -
  4718. -/* iv_add -- integer vector addition -- may be in-situ */
  4719. -IVEC    *iv_add(iv1,iv2,out)
  4720. -IVEC    *iv1,*iv2,*out;
  4721. -{
  4722. -   u_int    i;
  4723. -   int    *out_ive, *iv1_ive, *iv2_ive;
  4724. -   
  4725. -   if ( iv1==IVNULL || iv2==IVNULL )
  4726. -     error(E_NULL,"iv_add");
  4727. -   if ( iv1->dim != iv2->dim )
  4728. -     error(E_SIZES,"iv_add");
  4729. -   if ( out==IVNULL || out->dim != iv1->dim )
  4730. -     out = iv_resize(out,iv1->dim);
  4731. -   
  4732. -   out_ive = out->ive;
  4733. -   iv1_ive = iv1->ive;
  4734. -   iv2_ive = iv2->ive;
  4735. -   
  4736. -   for ( i = 0; i < iv1->dim; i++ )
  4737. -     out_ive[i] = iv1_ive[i] + iv2_ive[i];
  4738. -   
  4739. -   return (out);
  4740. -}
  4741. -
  4742. -
  4743. -
  4744. -/* iv_sub -- integer vector addition -- may be in-situ */
  4745. -IVEC    *iv_sub(iv1,iv2,out)
  4746. -IVEC    *iv1,*iv2,*out;
  4747. -{
  4748. -   u_int    i;
  4749. -   int    *out_ive, *iv1_ive, *iv2_ive;
  4750. -   
  4751. -   if ( iv1==IVNULL || iv2==IVNULL )
  4752. -     error(E_NULL,"iv_sub");
  4753. -   if ( iv1->dim != iv2->dim )
  4754. -     error(E_SIZES,"iv_sub");
  4755. -   if ( out==IVNULL || out->dim != iv1->dim )
  4756. -     out = iv_resize(out,iv1->dim);
  4757. -   
  4758. -   out_ive = out->ive;
  4759. -   iv1_ive = iv1->ive;
  4760. -   iv2_ive = iv2->ive;
  4761. -   
  4762. -   for ( i = 0; i < iv1->dim; i++ )
  4763. -     out_ive[i] = iv1_ive[i] - iv2_ive[i];
  4764. -   
  4765. -   return (out);
  4766. -}
  4767. -
  4768. -/* iv_foutput -- print a representation of iv on stream fp */
  4769. -void    iv_foutput(fp,iv)
  4770. -FILE    *fp;
  4771. -IVEC    *iv;
  4772. -{
  4773. -   int    i;
  4774. -   
  4775. -   fprintf(fp,"IntVector: ");
  4776. -   if ( iv == IVNULL )
  4777. -   {
  4778. -      fprintf(fp,"**** NULL ****\n");
  4779. -      return;
  4780. -   }
  4781. -   fprintf(fp,"dim: %d\n",iv->dim);
  4782. -   for ( i = 0; i < iv->dim; i++ )
  4783. -   {
  4784. -      if ( (i+1) % 8 )
  4785. -    fprintf(fp,"%8d ",iv->ive[i]);
  4786. -      else
  4787. -    fprintf(fp,"%8d\n",iv->ive[i]);
  4788. -   }
  4789. -   if ( i % 8 )
  4790. -     fprintf(fp,"\n");
  4791. -}
  4792. -
  4793. -
  4794. -/* iv_finput -- input integer vector from stream fp */
  4795. -IVEC    *iv_finput(fp,x)
  4796. -FILE    *fp;
  4797. -IVEC    *x;
  4798. -{
  4799. -   IVEC    *iiv_finput(),*biv_finput();
  4800. -   
  4801. -   if ( isatty(fileno(fp)) )
  4802. -     return iiv_finput(fp,x);
  4803. -   else
  4804. -     return biv_finput(fp,x);
  4805. -}
  4806. -
  4807. -/* iiv_finput -- interactive input of IVEC iv */
  4808. -IVEC    *iiv_finput(fp,iv)
  4809. -FILE    *fp;
  4810. -IVEC    *iv;
  4811. -{
  4812. -   u_int    i,dim,dynamic;    /* dynamic set if memory allocated here */
  4813. -   
  4814. -   /* get dimension */
  4815. -   if ( iv != (IVEC *)NULL && iv->dim<MAXDIM )
  4816. -   {    dim = iv->dim;    dynamic = FALSE;    }
  4817. -   else
  4818. -   {
  4819. -      dynamic = TRUE;
  4820. -      do
  4821. -      {
  4822. -     fprintf(stderr,"IntVector: dim: ");
  4823. -     if ( fgets(line,MAXLINE,fp)==NULL )
  4824. -       error(E_INPUT,"iiv_finput");
  4825. -      } while ( sscanf(line,"%u",&dim)<1 || dim>MAXDIM );
  4826. -      iv = iv_get(dim);
  4827. -   }
  4828. -   
  4829. -   /* input elements */
  4830. -   for ( i=0; i<dim; i++ )
  4831. -     do
  4832. -     {
  4833. -      redo:
  4834. -    fprintf(stderr,"entry %u: ",i);
  4835. -    if ( !dynamic )
  4836. -      fprintf(stderr,"old: %-9d  new: ",iv->ive[i]);
  4837. -    if ( fgets(line,MAXLINE,fp)==NULL )
  4838. -      error(E_INPUT,"iiv_finput");
  4839. -    if ( (*line == 'b' || *line == 'B') && i > 0 )
  4840. -    {    i--;    dynamic = FALSE;    goto redo;       }
  4841. -    if ( (*line == 'f' || *line == 'F') && i < dim-1 )
  4842. -    {    i++;    dynamic = FALSE;    goto redo;       }
  4843. -     } while ( *line=='\0' || sscanf(line,"%d",&iv->ive[i]) < 1 );
  4844. -   
  4845. -   return (iv);
  4846. -}
  4847. -
  4848. -/* biv_finput -- batch-file input of IVEC iv */
  4849. -IVEC    *biv_finput(fp,iv)
  4850. -FILE    *fp;
  4851. -IVEC    *iv;
  4852. -{
  4853. -   u_int    i,dim;
  4854. -   int    io_code;
  4855. -   
  4856. -   /* get dimension */
  4857. -   skipjunk(fp);
  4858. -   if ((io_code=fscanf(fp," IntVector: dim:%u",&dim)) < 1 ||
  4859. -       dim>MAXDIM )
  4860. -     error(io_code==EOF ? 7 : 6,"biv_finput");
  4861. -   
  4862. -   /* allocate memory if necessary */
  4863. -   if ( iv==(IVEC *)NULL || iv->dim<dim )
  4864. -     iv = iv_resize(iv,dim);
  4865. -   
  4866. -   /* get entries */
  4867. -   skipjunk(fp);
  4868. -   for ( i=0; i<dim; i++ )
  4869. -     if ((io_code=fscanf(fp,"%d",&iv->ive[i])) < 1 )
  4870. -       error(io_code==EOF ? 7 : 6,"biv_finput");
  4871. -   
  4872. -   return (iv);
  4873. -}
  4874. -
  4875. -/* iv_dump -- dumps all the contents of IVEC iv onto stream fp */
  4876. -void    iv_dump(fp,iv)
  4877. -FILE*fp;
  4878. -IVEC*iv;
  4879. -{
  4880. -   int        i;
  4881. -   
  4882. -   fprintf(fp,"IntVector: ");
  4883. -   if ( ! iv )
  4884. -   {
  4885. -      fprintf(fp,"**** NULL ****\n");
  4886. -      return;
  4887. -   }
  4888. -   fprintf(fp,"dim: %d, max_dim: %d\n",iv->dim,iv->max_dim);
  4889. -   fprintf(fp,"ive @ 0x%lx\n",(long)(iv->ive));
  4890. -   for ( i = 0; i < iv->max_dim; i++ )
  4891. -   {
  4892. -      if ( (i+1) % 8 )
  4893. -    fprintf(fp,"%8d ",iv->ive[i]);
  4894. -      else
  4895. -    fprintf(fp,"%8d\n",iv->ive[i]);
  4896. -   }
  4897. -   if ( i % 8 )
  4898. -     fprintf(fp,"\n");
  4899. -}    
  4900. -
  4901. -#define    MAX_STACK    60
  4902. -
  4903. -
  4904. -/* iv_sort -- sorts vector x, and generates permutation that gives the order
  4905. -   of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
  4906. -   the permutation is order = [2, 0, 1].
  4907. -   -- if order is NULL on entry then it is ignored
  4908. -   -- the sorted vector x is returned */
  4909. -IVEC    *iv_sort(x, order)
  4910. -IVEC    *x;
  4911. -PERM    *order;
  4912. -{
  4913. -   int        *x_ive, tmp, v;
  4914. -   /* int        *order_pe; */
  4915. -   int        dim, i, j, l, r, tmp_i;
  4916. -   int        stack[MAX_STACK], sp;
  4917. -   
  4918. -   if ( ! x )
  4919. -     error(E_NULL,"v_sort");
  4920. -   if ( order != PNULL && order->size != x->dim )
  4921. -     order = px_resize(order, x->dim);
  4922. -   
  4923. -   x_ive = x->ive;
  4924. -   dim = x->dim;
  4925. -   if ( order != PNULL )
  4926. -     px_ident(order);
  4927. -   
  4928. -   if ( dim <= 1 )
  4929. -     return x;
  4930. -   
  4931. -   /* using quicksort algorithm in Sedgewick,
  4932. -      "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
  4933. -   sp = 0;
  4934. -   l = 0;    r = dim-1;    v = x_ive[0];
  4935. -   for ( ; ; )
  4936. -   {
  4937. -      while ( r > l )
  4938. -      {
  4939. -     /* "i = partition(x_ive,l,r);" */
  4940. -     v = x_ive[r];
  4941. -     i = l-1;
  4942. -     j = r;
  4943. -     for ( ; ; )
  4944. -     {
  4945. -        while ( x_ive[++i] < v )
  4946. -          ;
  4947. -        while ( x_ive[--j] > v )
  4948. -          ;
  4949. -        if ( i >= j )    break;
  4950. -        
  4951. -        tmp = x_ive[i];
  4952. -        x_ive[i] = x_ive[j];
  4953. -        x_ive[j] = tmp;
  4954. -        if ( order != PNULL )
  4955. -        {
  4956. -           tmp_i = order->pe[i];
  4957. -           order->pe[i] = order->pe[j];
  4958. -           order->pe[j] = tmp_i;
  4959. -        }
  4960. -     }
  4961. -     tmp = x_ive[i];
  4962. -     x_ive[i] = x_ive[r];
  4963. -     x_ive[r] = tmp;
  4964. -     if ( order != PNULL )
  4965. -     {
  4966. -        tmp_i = order->pe[i];
  4967. -        order->pe[i] = order->pe[r];
  4968. -        order->pe[r] = tmp_i;
  4969. -     }
  4970. -     
  4971. -     if ( i-l > r-i )
  4972. -     {   stack[sp++] = l;   stack[sp++] = i-1;   l = i+1;   }
  4973. -     else
  4974. -     {   stack[sp++] = i+1;   stack[sp++] = r;   r = i-1;   }
  4975. -      }
  4976. -      
  4977. -      /* recursion elimination */
  4978. -      if ( sp == 0 )
  4979. -    break;
  4980. -      r = stack[--sp];
  4981. -      l = stack[--sp];
  4982. -   }
  4983. -   
  4984. -   return x;
  4985. -}
  4986. //GO.SYSIN DD ivecop.c
  4987. echo version.c 1>&2
  4988. sed >version.c <<'//GO.SYSIN DD version.c' 's/^-//'
  4989. -
  4990. -/**************************************************************************
  4991. -**
  4992. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4993. -**
  4994. -**                 Meschach Library
  4995. -** 
  4996. -** This Meschach Library is provided "as is" without any express 
  4997. -** or implied warranty of any kind with respect to this software. 
  4998. -** In particular the authors shall not be liable for any direct, 
  4999. -** indirect, special, incidental or consequential damages arising 
  5000. -** in any way from use of the software.
  5001. -** 
  5002. -** Everyone is granted permission to copy, modify and redistribute this
  5003. -** Meschach Library, provided:
  5004. -**  1.  All copies contain this copyright notice.
  5005. -**  2.  All modified copies shall carry a notice stating who
  5006. -**      made the last modification and the date of such modification.
  5007. -**  3.  No charge is made for this software or works derived from it.  
  5008. -**      This clause shall not be construed as constraining other software
  5009. -**      distributed on the same medium as this software, nor is a
  5010. -**      distribution fee considered a charge.
  5011. -**
  5012. -***************************************************************************/
  5013. -
  5014. -
  5015. -/*            Version routine            */
  5016. -/*    This routine must be modified whenever modifications are made to
  5017. -    Meschach by persons other than the original authors
  5018. -    (David E. Stewart & Zbigniew Leyk); 
  5019. -    when new releases of Meschach are made the
  5020. -    version number will also be updated
  5021. -*/
  5022. -
  5023. -#include    <stdio.h>
  5024. -
  5025. -void    m_version()
  5026. -{
  5027. -    static char rcsid[] = "$Id: version.c,v 1.9 1994/03/24 00:04:05 des Exp $";
  5028. -
  5029. -    printf("Meshach matrix library version 1.2b\n");
  5030. -    printf("RCS id: %s\n",rcsid);
  5031. -    printf("Changes since 1.2a:\n");
  5032. -    printf("\t Fixed bug in schur() for 2x2 blocks with real e-vals\n");
  5033. -    printf("\t Fixed bug in schur() reading beyond end of array\n");
  5034. -    printf("\t Fixed some installation bugs\n");
  5035. -    printf("\t Fixed bugs & improved efficiency in spILUfactor()\n");
  5036. -    printf("\t px_inv() doesn't crash inverting non-permutations\n");
  5037. -    /**** List of modifications ****/
  5038. -    /* Example below is for illustration only */
  5039. -    /* printf("Modified by %s, routine(s) %s, file %s on date %s\n",
  5040. -            "Joe Bloggs",
  5041. -            "m_version",
  5042. -            "version.c",
  5043. -            "Fri Apr  5 16:00:38 EST 1994"); */
  5044. -    /* printf("Purpose: %s\n",
  5045. -            "To update the version number"); */
  5046. -}
  5047. -
  5048. -/* $Log: version.c,v $
  5049. - * Revision 1.9  1994/03/24  00:04:05  des
  5050. - * Added notes on changes to spILUfactor() and px_inv().
  5051. - *
  5052. - * Revision 1.8  1994/02/21  04:32:25  des
  5053. - * Set version to 1.2b with bug fixes in schur() and installation.
  5054. - *
  5055. - * Revision 1.7  1994/01/13  05:43:57  des
  5056. - * Version 1.2 update
  5057. - *
  5058. -
  5059. - * */
  5060. //GO.SYSIN DD version.c
  5061. echo meminfo.c 1>&2
  5062. sed >meminfo.c <<'//GO.SYSIN DD meminfo.c' 's/^-//'
  5063. -
  5064. -/**************************************************************************
  5065. -**
  5066. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5067. -**
  5068. -**                 Meschach Library
  5069. -** 
  5070. -** This Meschach Library is provided "as is" without any express 
  5071. -** or implied warranty of any kind with respect to this software. 
  5072. -** In particular the authors shall not be liable for any direct, 
  5073. -** indirect, special, incidental or consequential damages arising 
  5074. -** in any way from use of the software.
  5075. -** 
  5076. -** Everyone is granted permission to copy, modify and redistribute this
  5077. -** Meschach Library, provided:
  5078. -**  1.  All copies contain this copyright notice.
  5079. -**  2.  All modified copies shall carry a notice stating who
  5080. -**      made the last modification and the date of such modification.
  5081. -**  3.  No charge is made for this software or works derived from it.  
  5082. -**      This clause shall not be construed as constraining other software
  5083. -**      distributed on the same medium as this software, nor is a
  5084. -**      distribution fee considered a charge.
  5085. -**
  5086. -***************************************************************************/
  5087. -
  5088. -
  5089. -/* meminfo.c  revised  22/11/93 */
  5090. -
  5091. -/* 
  5092. -  contains basic functions, types and arrays 
  5093. -  to keep track of memory allocation/deallocation
  5094. -*/
  5095. -
  5096. -#include <stdio.h>
  5097. -#include  "matrix.h"
  5098. -#include  "meminfo.h"
  5099. -#ifdef COMPLEX   
  5100. -#include  "zmatrix.h"
  5101. -#endif
  5102. -#ifdef SPARSE
  5103. -#include  "sparse.h"
  5104. -#include  "iter.h"
  5105. -#endif
  5106. -
  5107. -static char rcsid[] = "$Id: meminfo.c,v 1.1 1994/01/13 05:31:39 des Exp $";
  5108. -
  5109. -/* this array is defined further in this file */
  5110. -extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS];
  5111. -
  5112. -
  5113. -/* names of types */
  5114. -static char *mem_type_names[] = {
  5115. -   "MAT",
  5116. -   "BAND",
  5117. -   "PERM",
  5118. -   "VEC",
  5119. -   "IVEC"
  5120. -#ifdef SPARSE
  5121. -     ,"ITER",
  5122. -     "SPROW",
  5123. -     "SPMAT"
  5124. -#endif
  5125. -#ifdef COMPLEX   
  5126. -       ,"ZVEC",
  5127. -       "ZMAT"
  5128. -#endif
  5129. -      };
  5130. -
  5131. -
  5132. -#define MEM_NUM_STD_TYPES  (sizeof(mem_type_names)/sizeof(mem_type_names[0]))
  5133. -
  5134. -
  5135. -/* local array for keeping track of memory */
  5136. -static MEM_ARRAY   mem_info_sum[MEM_NUM_STD_TYPES];  
  5137. -
  5138. -
  5139. -/* for freeing various types */
  5140. -static int (*mem_free_funcs[MEM_NUM_STD_TYPES])() = {
  5141. -   m_free,
  5142. -   bd_free,
  5143. -   px_free,    
  5144. -   v_free,    
  5145. -   iv_free
  5146. -#ifdef SPARSE
  5147. -     ,iter_free,    
  5148. -     sprow_free, 
  5149. -     sp_free
  5150. -#endif
  5151. -#ifdef COMPLEX
  5152. -       ,zv_free,    
  5153. -       zm_free
  5154. -#endif
  5155. -      };
  5156. -
  5157. -
  5158. -
  5159. -/* it is a global variable for passing 
  5160. -   pointers to local arrays defined here */
  5161. -MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS] = {
  5162. - { mem_type_names, mem_free_funcs, MEM_NUM_STD_TYPES, 
  5163. -     mem_info_sum } 
  5164. -};
  5165. -
  5166. -
  5167. -/* attach a new list of types */
  5168. -
  5169. -int mem_attach_list(list, ntypes, type_names, free_funcs, info_sum)
  5170. -int list,ntypes;         /* number of a list and number of types there */
  5171. -char *type_names[];      /* list of names of types */
  5172. -int (*free_funcs[])();   /* list of releasing functions */
  5173. -MEM_ARRAY info_sum[];    /* local table */
  5174. -{
  5175. -   if (list < 0 || list >= MEM_CONNECT_MAX_LISTS)
  5176. -     return -1;
  5177. -
  5178. -   if (type_names == NULL || free_funcs == NULL 
  5179. -       || info_sum == NULL || ntypes < 0)
  5180. -     return -1;
  5181. -   
  5182. -   /* if a list exists do not overwrite */
  5183. -   if ( mem_connect[list].ntypes != 0 )
  5184. -     error(E_OVERWRITE,"mem_attach_list");
  5185. -   
  5186. -   mem_connect[list].ntypes = ntypes;
  5187. -   mem_connect[list].type_names = type_names;
  5188. -   mem_connect[list].free_funcs = free_funcs;
  5189. -   mem_connect[list].info_sum = info_sum;
  5190. -   return 0;
  5191. -}
  5192. -
  5193. -
  5194. -/* release a list of types */
  5195. -int mem_free_vars(list)
  5196. -int list;
  5197. -{    
  5198. -   if (list < 0 || list >= MEM_CONNECT_MAX_LISTS)
  5199. -     return -1;
  5200. -   
  5201. -   mem_connect[list].ntypes = 0;
  5202. -   mem_connect[list].type_names = NULL;
  5203. -   mem_connect[list].free_funcs = NULL;
  5204. -   mem_connect[list].info_sum = NULL;
  5205. -   
  5206. -   return 0;
  5207. -}
  5208. -
  5209. -
  5210. -
  5211. -/* check if list is attached */
  5212. -
  5213. -int mem_is_list_attached(list)
  5214. -int list;
  5215. -{
  5216. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5217. -   return FALSE;
  5218. -
  5219. -   if ( mem_connect[list].type_names != NULL &&
  5220. -        mem_connect[list].free_funcs != NULL &&
  5221. -        mem_connect[list].info_sum != NULL)
  5222. -     return TRUE;
  5223. -   else return FALSE;
  5224. -}
  5225. -
  5226. -/* to print out the contents of mem_connect[list] */
  5227. -
  5228. -void mem_dump_list(fp,list)
  5229. -FILE *fp;
  5230. -int list;
  5231. -{
  5232. -   int i;
  5233. -   MEM_CONNECT *mlist;
  5234. -
  5235. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5236. -     return;
  5237. -
  5238. -   mlist = &mem_connect[list];
  5239. -   fprintf(fp," %15s[%d]:\n","CONTENTS OF mem_connect",list);
  5240. -   fprintf(fp," %-7s   %-12s   %-9s   %s\n",
  5241. -       "name of",
  5242. -       "alloc.", "# alloc.",
  5243. -       "address"
  5244. -       );
  5245. -   fprintf(fp," %-7s   %-12s   %-9s   %s\n",
  5246. -       " type",
  5247. -       "bytes", "variables",
  5248. -       "of *_free()"
  5249. -       );
  5250. -
  5251. -   for (i=0; i < mlist->ntypes; i++) 
  5252. -     fprintf(fp,"  %-7s   %-12ld   %-9d   %p\n",
  5253. -         mlist->type_names[i], mlist->info_sum[i].bytes,
  5254. -         mlist->info_sum[i].numvar, mlist->free_funcs[i]
  5255. -         );
  5256. -   
  5257. -   fprintf(fp,"\n");
  5258. -}
  5259. -
  5260. -
  5261. -
  5262. -/*=============================================================*/
  5263. -
  5264. -
  5265. -/* local variables */
  5266. -
  5267. -static int    mem_switched_on = MEM_SWITCH_ON_DEF;  /* on/off */
  5268. -
  5269. -
  5270. -/* switch on/off memory info */
  5271. -
  5272. -int mem_info_on(sw)
  5273. -int sw;
  5274. -{
  5275. -   int old = mem_switched_on;
  5276. -   
  5277. -   mem_switched_on = sw;
  5278. -   return old;
  5279. -}
  5280. -
  5281. -#ifdef ANSI_C
  5282. -int mem_info_is_on(void)
  5283. -#else
  5284. -int mem_info_is_on()
  5285. -#endif
  5286. -{
  5287. -   return mem_switched_on;
  5288. -}
  5289. -
  5290. -
  5291. -/* information about allocated memory */
  5292. -
  5293. -/* return the number of allocated bytes for type 'type' */
  5294. -
  5295. -long mem_info_bytes(type,list)
  5296. -int type,list;
  5297. -{
  5298. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5299. -     return 0l;
  5300. -   if ( !mem_switched_on || type < 0 
  5301. -       || type >= mem_connect[list].ntypes
  5302. -       || mem_connect[list].free_funcs[type] == NULL )
  5303. -     return 0l;
  5304. -   
  5305. -   return mem_connect[list].info_sum[type].bytes;
  5306. -}
  5307. -
  5308. -/* return the number of allocated variables for type 'type' */
  5309. -int mem_info_numvar(type,list)
  5310. -int type,list;
  5311. -{
  5312. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5313. -     return 0l;
  5314. -   if ( !mem_switched_on || type < 0 
  5315. -       || type >= mem_connect[list].ntypes
  5316. -       || mem_connect[list].free_funcs[type] == NULL )
  5317. -     return 0l;
  5318. -   
  5319. -   return mem_connect[list].info_sum[type].numvar;
  5320. -}
  5321. -
  5322. -
  5323. -
  5324. -/* print out memory info to the file fp */
  5325. -void mem_info_file(fp,list)
  5326. -FILE *fp;
  5327. -int list;
  5328. -{
  5329. -   unsigned int type;
  5330. -   long t = 0l, d;
  5331. -   int n = 0, nt = 0;
  5332. -   MEM_CONNECT *mlist;
  5333. -   
  5334. -   if (!mem_switched_on) return;
  5335. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5336. -     return;
  5337. -   
  5338. -   if (list == 0)
  5339. -     fprintf(fp," MEMORY INFORMATION (standard types):\n");
  5340. -   else
  5341. -     fprintf(fp," MEMORY INFORMATION (list no. %d):\n",list);
  5342. -
  5343. -   mlist = &mem_connect[list];
  5344. -
  5345. -   for (type=0; type < mlist->ntypes; type++) {
  5346. -      if (mlist->type_names[type] == NULL ) continue;
  5347. -      d = mlist->info_sum[type].bytes;
  5348. -      t += d;
  5349. -      n = mlist->info_sum[type].numvar;
  5350. -      nt += n;
  5351. -      fprintf(fp," type %-7s %10ld alloc. byte%c  %6d alloc. variable%c\n",
  5352. -          mlist->type_names[type], d, (d!=1 ? 's' : ' '),
  5353. -          n, (n!=1 ? 's' : ' '));
  5354. -   }
  5355. -
  5356. -   fprintf(fp," %-12s %10ld alloc. byte%c  %6d alloc. variable%c\n\n",
  5357. -       "total:",t, (t!=1 ? 's' : ' '),
  5358. -       nt, (nt!=1 ? 's' : ' '));
  5359. -}
  5360. -
  5361. -
  5362. -/* function for memory information */
  5363. -
  5364. -
  5365. -/* mem_bytes_list
  5366. -   
  5367. -   Arguments:
  5368. -   type - the number of type;
  5369. -   old_size - old size of allocated memory (in bytes);
  5370. -   new_size - new size of allocated memory (in bytes);
  5371. -   list - list of types
  5372. -   */
  5373. -
  5374. -
  5375. -void mem_bytes_list(type,old_size,new_size,list)
  5376. -int type,list;
  5377. -int old_size,new_size;
  5378. -{
  5379. -   MEM_CONNECT *mlist;
  5380. -   
  5381. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5382. -     return;
  5383. -   
  5384. -   mlist = &mem_connect[list];
  5385. -   if (  type < 0 || type >= mlist->ntypes
  5386. -       || mlist->free_funcs[type] == NULL )
  5387. -     return;
  5388. -
  5389. -   if ( old_size < 0 || new_size < 0 )
  5390. -     error(E_NEG,"mem_bytes_list");
  5391. -
  5392. -   mlist->info_sum[type].bytes += new_size - old_size;
  5393. -   
  5394. -   /* check if the number of bytes is non-negative */
  5395. -   if ( old_size > 0 ) {
  5396. -
  5397. -      if (mlist->info_sum[type].bytes < 0)
  5398. -      {
  5399. -     fprintf(stderr,
  5400. -       "\n WARNING !! memory info: allocated memory is less than 0\n");
  5401. -     fprintf(stderr,"\t TYPE %s \n\n", mlist->type_names[type]);
  5402. -
  5403. -     if ( !isatty(fileno(stdout)) ) {
  5404. -        fprintf(stdout,
  5405. -          "\n WARNING !! memory info: allocated memory is less than 0\n");
  5406. -        fprintf(stdout,"\t TYPE %s \n\n", mlist->type_names[type]);
  5407. -     }
  5408. -      }
  5409. -   }
  5410. -}
  5411. -
  5412. -
  5413. -/* mem_numvar_list
  5414. -   
  5415. -   Arguments:
  5416. -   type - the number of type;
  5417. -   num - # of variables allocated (> 0) or deallocated ( < 0)
  5418. -   list - list of types
  5419. -   */
  5420. -
  5421. -
  5422. -void mem_numvar_list(type,num,list)
  5423. -int type,list,num;
  5424. -{
  5425. -   MEM_CONNECT *mlist;
  5426. -   
  5427. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5428. -     return;
  5429. -   
  5430. -   mlist = &mem_connect[list];
  5431. -   if (  type < 0 || type >= mlist->ntypes
  5432. -       || mlist->free_funcs[type] == NULL )
  5433. -     return;
  5434. -
  5435. -   mlist->info_sum[type].numvar += num;
  5436. -   
  5437. -   /* check if the number of variables is non-negative */
  5438. -   if ( num < 0 ) {
  5439. -
  5440. -      if (mlist->info_sum[type].numvar < 0)
  5441. -      {
  5442. -     fprintf(stderr,
  5443. -       "\n WARNING !! memory info: allocated # of variables is less than 0\n");
  5444. -     fprintf(stderr,"\t TYPE %s \n\n", mlist->type_names[type]);
  5445. -     if ( !isatty(fileno(stdout)) ) {
  5446. -        fprintf(stdout,
  5447. -      "\n WARNING !! memory info: allocated # of variables is less than 0\n");
  5448. -        fprintf(stdout,"\t TYPE %s \n\n", mlist->type_names[type]);
  5449. -     }
  5450. -      }
  5451. -   }
  5452. -}
  5453. -
  5454. //GO.SYSIN DD meminfo.c
  5455. echo memstat.c 1>&2
  5456. sed >memstat.c <<'//GO.SYSIN DD memstat.c' 's/^-//'
  5457. -
  5458. -/**************************************************************************
  5459. -**
  5460. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5461. -**
  5462. -**                 Meschach Library
  5463. -** 
  5464. -** This Meschach Library is provided "as is" without any express 
  5465. -** or implied warranty of any kind with respect to this software. 
  5466. -** In particular the authors shall not be liable for any direct, 
  5467. -** indirect, special, incidental or consequential damages arising 
  5468. -** in any way from use of the software.
  5469. -** 
  5470. -** Everyone is granted permission to copy, modify and redistribute this
  5471. -** Meschach Library, provided:
  5472. -**  1.  All copies contain this copyright notice.
  5473. -**  2.  All modified copies shall carry a notice stating who
  5474. -**      made the last modification and the date of such modification.
  5475. -**  3.  No charge is made for this software or works derived from it.  
  5476. -**      This clause shall not be construed as constraining other software
  5477. -**      distributed on the same medium as this software, nor is a
  5478. -**      distribution fee considered a charge.
  5479. -**
  5480. -***************************************************************************/
  5481. -
  5482. -
  5483. -/*  mem_stat.c    6/09/93  */
  5484. -
  5485. -/* Deallocation of static arrays */
  5486. -
  5487. -
  5488. -#include <stdio.h>
  5489. -#include  "matrix.h"
  5490. -#include  "meminfo.h"
  5491. -#ifdef COMPLEX   
  5492. -#include  "zmatrix.h"
  5493. -#endif
  5494. -#ifdef SPARSE
  5495. -#include  "sparse.h"
  5496. -#include  "iter.h"
  5497. -#endif
  5498. -
  5499. -static char rcsid[] = "$Id: memstat.c,v 1.1 1994/01/13 05:32:44 des Exp $";
  5500. -
  5501. -/* global variable */
  5502. -
  5503. -extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS];
  5504. -
  5505. -
  5506. -/* local type */
  5507. -
  5508. -typedef struct {
  5509. -   void **var;    /* for &A, where A is a pointer */
  5510. -   int type;     /* type of A */
  5511. -   int mark;      /* what mark is chosen */
  5512. -} MEM_STAT_STRUCT;
  5513. -
  5514. -
  5515. -/* local variables */
  5516. -
  5517. -/* how many marks are used */
  5518. -static int mem_stat_mark_many = 0;
  5519. -
  5520. -/* current mark */
  5521. -static int mem_stat_mark_curr = 0;
  5522. -
  5523. -
  5524. -static MEM_STAT_STRUCT mem_stat_var[MEM_HASHSIZE];
  5525. -
  5526. -/* array of indices (+1) to mem_stat_var */
  5527. -static unsigned int mem_hash_idx[MEM_HASHSIZE];
  5528. -
  5529. -/* points to the first unused element in mem_hash_idx */
  5530. -static unsigned int mem_hash_idx_end = 0;
  5531. -
  5532. -
  5533. -
  5534. -/* hashing function */
  5535. -
  5536. -static unsigned int mem_hash(ptr)
  5537. -void **ptr;
  5538. -{
  5539. -   unsigned long lp = (unsigned long)ptr;
  5540. -
  5541. -   return (lp % MEM_HASHSIZE);
  5542. -}
  5543. -
  5544. -
  5545. -/* look for a place in mem_stat_var */
  5546. -static int mem_lookup(var)
  5547. -void **var;
  5548. -{
  5549. -   int k, j;
  5550. -
  5551. -   k = mem_hash(var);
  5552. -
  5553. -   if (mem_stat_var[k].var == var) {
  5554. -      return -1;
  5555. -   }
  5556. -   else if (mem_stat_var[k].var == NULL) {
  5557. -      return k;
  5558. -   }
  5559. -   else {  /* look for an empty place */
  5560. -      j = k;
  5561. -      while (mem_stat_var[j].var != var && j < MEM_HASHSIZE
  5562. -         && mem_stat_var[j].var != NULL) 
  5563. -    j++;
  5564. -
  5565. -      if (mem_stat_var[j].var == NULL) return j;
  5566. -      else if (mem_stat_var[j].var == var) return -1; 
  5567. -      else { /* if (j == MEM_HASHSIZE) */
  5568. -     j = 0;
  5569. -     while (mem_stat_var[j].var != var && j < k
  5570. -        && mem_stat_var[j].var != NULL) 
  5571. -       j++;
  5572. -     if (mem_stat_var[j].var == NULL) return j;
  5573. -     else if (mem_stat_var[j].var == var) return -1; 
  5574. -     else { /* if (j == k) */
  5575. -        fprintf(stderr,
  5576. -              "\n WARNING !!! static memory: mem_stat_var is too small\n");
  5577. -        fprintf(stderr,
  5578. -          " Increase MEM_HASHSIZE in file: %s (currently = %d)\n\n",
  5579. -            MEM_HASHSIZE_FILE, MEM_HASHSIZE);
  5580. -        if ( !isatty(fileno(stdout)) ) {
  5581. -           fprintf(stdout,
  5582. -                "\n WARNING !!! static memory: mem_stat_var is too small\n");
  5583. -           fprintf(stdout,
  5584. -            " Increase MEM_HASHSIZE in file: %s (currently = %d)\n\n",
  5585. -            MEM_HASHSIZE_FILE, MEM_HASHSIZE);
  5586. -        }
  5587. -        error(E_MEM,"mem_lookup");
  5588. -     }
  5589. -      }
  5590. -   }
  5591. -
  5592. -   return -1;
  5593. -}
  5594. -
  5595. -
  5596. -/* register static variables;
  5597. -   Input arguments:
  5598. -     var - variable to be registered,
  5599. -     type - type of this variable; 
  5600. -     list - list of types
  5601. -
  5602. -   returned value < 0  --> error,
  5603. -   returned value == 0 --> not registered,
  5604. -   returned value >= 0 --> registered with this mark;
  5605. -*/
  5606. -
  5607. -int mem_stat_reg_list(var,type,list)
  5608. -void **var;
  5609. -int type,list;
  5610. -{
  5611. -   int n;
  5612. -
  5613. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS )
  5614. -     return -1;
  5615. -
  5616. -   if (mem_stat_mark_curr == 0) return 0;  /* not registered */
  5617. -   if (var == NULL) return -1;             /* error */
  5618. -
  5619. -   if ( type < 0 || type >= mem_connect[list].ntypes || 
  5620. -       mem_connect[list].free_funcs[type] == NULL )
  5621. -   {
  5622. -      warning(WARN_WRONG_TYPE,"mem_stat_reg_list");
  5623. -      return -1;
  5624. -   }
  5625. -   
  5626. -   if ((n = mem_lookup(var)) >= 0) {
  5627. -      mem_stat_var[n].var = var;
  5628. -      mem_stat_var[n].mark = mem_stat_mark_curr;
  5629. -      mem_stat_var[n].type = type;
  5630. -      /* save n+1, not n */
  5631. -      mem_hash_idx[mem_hash_idx_end++] = n+1;
  5632. -   }
  5633. -
  5634. -   return mem_stat_mark_curr;
  5635. -}
  5636. -
  5637. -
  5638. -/* set a mark;
  5639. -   Input argument:
  5640. -   mark - positive number denoting a mark;
  5641. -   returned: 
  5642. -             mark if mark > 0,
  5643. -             0 if mark == 0,
  5644. -         -1 if mark is negative.
  5645. -*/
  5646. -
  5647. -int mem_stat_mark(mark)
  5648. -int mark;
  5649. -{
  5650. -   if (mark < 0) {
  5651. -      mem_stat_mark_curr = 0;
  5652. -      return -1;   /* error */
  5653. -   }
  5654. -   else if (mark == 0) {
  5655. -      mem_stat_mark_curr = 0; 
  5656. -      return 0; 
  5657. -   }
  5658. -
  5659. -   mem_stat_mark_curr = mark;
  5660. -   mem_stat_mark_many++;
  5661. -
  5662. -   return mark;
  5663. -}
  5664. -
  5665. -
  5666. -
  5667. -/* deallocate static variables;
  5668. -   Input argument:
  5669. -   mark - a positive number denoting the mark;
  5670. -
  5671. -   Returned:
  5672. -     -1 if mark < 0 (error);
  5673. -     0  if mark == 0;
  5674. -*/
  5675. -
  5676. -int mem_stat_free_list(mark,list)
  5677. -int mark,list;
  5678. -{
  5679. -   u_int i,j;
  5680. -   int     (*free_fn)();
  5681. -
  5682. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS 
  5683. -       || mem_connect[list].free_funcs == NULL )
  5684. -     return -1;
  5685. -
  5686. -   if (mark < 0) {
  5687. -      mem_stat_mark_curr = 0;
  5688. -      return -1;
  5689. -   }
  5690. -   else if (mark == 0) {
  5691. -      mem_stat_mark_curr = 0;
  5692. -      return 0;
  5693. -   }
  5694. -   
  5695. -   if (mem_stat_mark_many <= 0) {
  5696. -      warning(WARN_NO_MARK,"mem_stat_free");
  5697. -      return -1;
  5698. -   }
  5699. -
  5700. -   /* deallocate the marked variables */
  5701. -   for (i=0; i < mem_hash_idx_end; i++) {
  5702. -      j = mem_hash_idx[i];
  5703. -      if (j == 0) continue;
  5704. -      else {
  5705. -     j--;
  5706. -     if (mem_stat_var[j].mark == mark) {
  5707. -         free_fn = mem_connect[list].free_funcs[mem_stat_var[j].type];
  5708. -         if ( free_fn != NULL )
  5709. -         (*free_fn)(*mem_stat_var[j].var);
  5710. -         else
  5711. -         warning(WARN_WRONG_TYPE,"mem_stat_free");
  5712. -        
  5713. -        *(mem_stat_var[j].var) = NULL;
  5714. -        mem_stat_var[j].var = NULL;
  5715. -        mem_stat_var[j].mark = 0;
  5716. -        mem_hash_idx[i] = 0;
  5717. -     }
  5718. -      }
  5719. -   }
  5720. -
  5721. -   while (mem_hash_idx_end > 0 && mem_hash_idx[mem_hash_idx_end-1] == 0)
  5722. -     mem_hash_idx_end--;
  5723. -
  5724. -   mem_stat_mark_curr = 0;
  5725. -   mem_stat_mark_many--;
  5726. -   return 0;
  5727. -}
  5728. -
  5729. -
  5730. -/* only for diagnostic purposes */
  5731. -
  5732. -void mem_stat_dump(fp,list)
  5733. -FILE *fp;
  5734. -int list;
  5735. -{
  5736. -   u_int i,j,k=1;
  5737. -
  5738. -   if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS 
  5739. -       || mem_connect[list].free_funcs == NULL )
  5740. -     return;
  5741. -
  5742. -   fprintf(fp," Array mem_stat_var (list no. %d):\n",list);
  5743. -   for (i=0; i < mem_hash_idx_end; i++) {
  5744. -      j = mem_hash_idx[i];
  5745. -      if (j == 0) continue;
  5746. -      else {
  5747. -     j--;
  5748. -     fprintf(fp," %d.  var = 0x%p, type = %s, mark = %d\n",
  5749. -         k,mem_stat_var[j].var,
  5750. -         mem_stat_var[j].type < mem_connect[list].ntypes &&
  5751. -         mem_connect[list].free_funcs[mem_stat_var[j].type] != NULL ?
  5752. -         mem_connect[list].type_names[(int)mem_stat_var[j].type] : 
  5753. -         "???",
  5754. -         mem_stat_var[j].mark);
  5755. -     k++;
  5756. -      }
  5757. -   }
  5758. -   
  5759. -   fprintf(fp,"\n");
  5760. -}
  5761. -
  5762. -
  5763. -/* query function about the current mark */
  5764. -#ifdef ANSI_C
  5765. -int mem_stat_show_mark(void)
  5766. -#else
  5767. -int mem_stat_show_mark()
  5768. -#endif
  5769. -{
  5770. -   return mem_stat_mark_curr;
  5771. -}
  5772. -
  5773. -
  5774. -/* Varying number of arguments */
  5775. -
  5776. -
  5777. -#ifdef ANSI_C
  5778. -
  5779. -/* To allocate memory to many arguments. 
  5780. -   The function should be called:
  5781. -   mem_stat_vars(list,type,&v1,&v2,&v3,...,VNULL);
  5782. -   where 
  5783. -     int list,type;
  5784. -     void **v1, **v2, **v3,...;
  5785. -     The last argument should be VNULL ! 
  5786. -     type is the type of variables v1,v2,v3,...
  5787. -     (of course they must be of the same type)
  5788. -*/
  5789. -
  5790. -int mem_stat_reg_vars(int list,int type,...)
  5791. -{
  5792. -   va_list ap;
  5793. -   int i=0;
  5794. -   void **par;
  5795. -   
  5796. -   va_start(ap, type);
  5797. -   while (par = va_arg(ap,void **)) {   /* NULL ends the list*/
  5798. -      mem_stat_reg_list(par,type,list);
  5799. -      i++;
  5800. -   } 
  5801. -
  5802. -   va_end(ap);
  5803. -   return i;
  5804. -}
  5805. -
  5806. -#elif VARARGS
  5807. -/* old varargs is used */
  5808. -
  5809. -/* To allocate memory to many arguments. 
  5810. -   The function should be called:
  5811. -   mem_stat_vars(list,type,&v1,&v2,&v3,...,VNULL);
  5812. -   where 
  5813. -     int list,type;
  5814. -     void **v1, **v2, **v3,...;
  5815. -     The last argument should be VNULL ! 
  5816. -     type is the type of variables v1,v2,v3,...
  5817. -     (of course they must be of the same type)
  5818. -*/
  5819. -
  5820. -int mem_stat_reg_vars(va_alist) va_dcl
  5821. -{
  5822. -   va_list ap;
  5823. -   int type,list,i=0;
  5824. -   void **par;
  5825. -   
  5826. -   va_start(ap);
  5827. -   list = va_arg(ap,int);
  5828. -   type = va_arg(ap,int);
  5829. -   while (par = va_arg(ap,void **)) {   /* NULL ends the list*/
  5830. -      mem_stat_reg_list(par,type,list);
  5831. -      i++;
  5832. -   } 
  5833. -
  5834. -   va_end(ap);
  5835. -   return i;
  5836. -}
  5837. -
  5838. -
  5839. -#endif
  5840. //GO.SYSIN DD memstat.c
  5841.